adpak.f

      SUBROUTINE ADCE
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     CALCULATE CORONAL EQUILIBRIUM (INCLUDING CHARGE
C     EXCHANGE RECOMBINATION, IF PRESENT) USING THE IONIZATION,
C     RECOMBINATION, AND RADIATION RATES
C     STORED IN THE APPROPRIATE ARRAYS BY PREVIOUS CALLS
C     TO ADSET, ADERC, AND ADBCXR.
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
      COMMON / ADNEUT / RCXREC(100), RADBCX(100)
C
      COMMON / ADCEQ / CEFRAC(100), CEAVGZ, CERAD
C
      DIMENSION ZIONT(100), ZRECT(100), ZRADT(100)
      DIMENSION ZI(1), ZRE(1), ZRA(1)
C
C     FILL TOTAL IONIZATION, RECOMBINATION, AND RADIATION ARRAYS
C
      DO 10 JSPC = 1, NSPC
C
      ZIONT(JSPC) = RCLION(JSPC)
      ZRECT(JSPC) = RRAREC(JSPC) + RDIREC(JSPC) + RCXREC(JSPC)
      ZRADT(JSPC) = RADRRC(JSPC) + RADDRC(JSPC) + RADCLX(JSPC)
     &                           + RADBRM(JSPC) + RADBCX(JSPC)
C
   10 CONTINUE
C
      GO TO 15
C
C     ENTRY FOR ALTERNATIVE EXTERNAL INPUT OF TOTAL RATE DATA
C
      ENTRY ADCE2(ZI, ZRE, ZRA)
C
      DO 12 JSPC= 1, NSPC
      ZIONT(JSPC) = ZI(JSPC)
      ZRECT(JSPC) = ZRE(JSPC)
      ZRADT(JSPC) = ZRA(JSPC)
   12 CONTINUE
C
C     FIRST GUESS AT MAXIMUM SPECIES
C
   15 DO 20 JQ = 1 , NSPC
      IF(ZRECT(JQ) .GE. ZIONT(JQ))GO TO 21
   20 CONTINUE
      JQ = NSPC
   21 ICNTR = JQ
CW    WRITE(IOUT,23)ICNTR
   23 FORMAT(1X,'ICNTR 1ST GUESS =',I3)
      IPASS = 1
C
C     INITIALIZE CEFRAC ARRAY
C
   25 DO 26 JQ = 1,NSPC
   26 CEFRAC(JQ) = 0.
      CEFRAC(ICNTR) = 1.
C
C     IF ICNTR IS FULLY STRIPPED WORK DOWN IN IONIZ ONLY
C
      IF(ICNTR.EQ.NSPC)GO TO 35
C
C     WORK UP FROM ICNTR SPECIES (TOWARDS HIGHER IONIZ STATES)
C
      II = ICNTR + 1
C
      DO 30 JQ = II,NSPC
      CEFRAC(JQ) = CEFRAC(JQ-1) * ZIONT(JQ-1) / ZRECT(JQ)
      IF(CEFRAC(JQ).LT.1.E-10*CEFRAC(ICNTR)) GO TO 32
   30 CONTINUE
C
C     IF ICNTR IS NEUTRAL WE'RE DONE
C
   32 IF(ICNTR.EQ.1)GO TO 42
C
C     WORK DOWN FROM ICNTR SPECIES
C
   35 II = ICNTR - 1
      DO 40 J = 1,II
      JQ = II - J + 1
      CEFRAC(JQ) = CEFRAC(JQ+1) * ZRECT(JQ+1) / ZIONT(JQ)
      IF(CEFRAC(JQ).LT.1.E-10*CEFRAC(ICNTR)) GO TO 42
   40 CONTINUE
C
C
   42 IF(IPASS.EQ.2)GO TO 49
C
C     END OF FIRST PASS - FIND NEW MAX SPECIES AND ITERATE
C
      IPASS = 2
      ZFRCMX = 0.
      DO 45 JQ = 1, NSPC
      IF(CEFRAC(JQ).LT.ZFRCMX)GO TO 45
      ZFRCMX = CEFRAC(JQ)
      IQFRAC = JQ
   45 CONTINUE
      ICNTR = IQFRAC
      GO TO 25
C
C
C
   49 ZPSUM = 0.
      DO 50 JQ = 1, NSPC
   50 ZPSUM = ZPSUM + CEFRAC(JQ)
C
C     NORMALIZE SPECIES FRACTIONS AND COMPUTE RADIATION RATE
C
C
      CERAD = 0.
C
      DO 60 JQ = 1, NSPC
      CEFRAC(JQ) = CEFRAC(JQ) / ZPSUM
      IF(ZRADT(JQ) .LE. 0.) GO TO 60
      IF(CEFRAC(JQ) .GT. 1.E-37 / ZRADT(JQ))
     &                       CERAD = CEFRAC(JQ) * ZRADT(JQ) + CERAD
C
   60 CONTINUE
C
C     COMPUTE <Z>
C
      CEAVGZ = 0.
      DO 70 JQ = 2,NSPC
   70 CEAVGZ = CEAVGZ + FLOAT(JQ-1) * CEFRAC(JQ)
C
C
      RETURN
      END
      SUBROUTINE ADSET(KNUCZ)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     SETUP BASIC ATOMIC STRUCTURE DATA FOR CHOSEN ELEMENT
C     (SHELL POPULATIONS, ENERGIES, ETC.)
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADPARM / LADEN, LADTIP, LTIPOK, LECI, LECIOK
C
      NUCZ = KNUCZ
C
      CALL ADAPN
C
      IF(LADEN .NE. 0 .AND. LADEN .NE. 1) STOP
      IF(LADEN .EQ. 0) CALL ADMAYR
      IF(LADEN .EQ. 1) CALL ADMORE
C
      IF(LADTIP .NE. 0 .AND. LADTIP .NE. 1) STOP
      LTIPOK = 0
      IF(LADTIP .EQ. 1) CALL ADTIP(LTIPOK)
C
      CALL ADFNM
C
      RETURN
      END
      SUBROUTINE ADAPN
C
C     REV: 10/9/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     INPUT : NUCZ
C     OUTPUT : NSPC, APN, AQN, NVALNC
C
C     ASSIGN THE NUMBER OF ELECTRONS PER SHELL FOR EACH
C     IONIC SPECIES.  APN(JQ,JN) = NUMBER OF ELECTRONS IN
C     SHELL JN FOR SPECIES WITH CHARGE JQ-1, WHILE AQN(JQ,JN) =
C     (2.*JN**2 - APN(JQ,JN))/(2.*JN**2) IS THE FRACTIONAL VACANCY.
C     THE SHELLS ARE FILLED IN SIMPLE CONSECUTIVE ASCENDING ORDER.
C     ALSO, THE PRINCIPAL QUANTUM NUMBER OF THE HIGHEST
C     OCCUPIED ('VALENCE') SHELL FOR EACH SPECIES IS STORED
C     IN NVALNC(JQ).
C
C     NOTE:  FOR THE FULLY IONIZED SPECIES NVALNC IS DEFAULTED TO 1.
C     ALSO: NSPC = NUCZ + 1 IS SET BY THIS ROUTINE
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
C     SET MAXIMUM SPECIES INDEX NSPC = NUCZ + 1
C
      NSPC = NUCZ +1
C
C
      DO 1 JQ = 1,NSPC
      ICOUNT = NUCZ-JQ+1
C
C
      DO 2 JN = 1,10
      IFULL = 2*JN*JN
      ISHELL = MIN0 (ICOUNT,IFULL)
      ICOUNT = ICOUNT - ISHELL
      IF(ISHELL.GT.0) NVALNC(JQ) = JN
      APN(JQ,JN) = FLOAT(ISHELL)
      AQN(JQ,JN) = FLOAT(IFULL - ISHELL)/FLOAT(IFULL)
C
CW    WRITE(6,900)JQ,JN,APN(JQ,JN),NVALNC(JQ)
CW  900 FORMAT(1X,'JQ =',I3,' JN =',I3,'  APN=',F4.0,'  NVALNC=',I3)
C
    2 CONTINUE
C
C
    1 CONTINUE
C
C
      NVALNC(NSPC) = 1
C
      RETURN
      END
      SUBROUTINE ADMORE
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 1 SEPTEMBER 1981
C
C     INPUT : NUCZ, NSPC, APN, NVALNC
C     OUTPUT : EN, EIN, ENM, QN
C
C     EXCITATION ENERGIES AND IONIZATION POTENTIALS ARE
C     TABULATED FOR ALL OCCUPIED SHELLS FOR EACH SPECIES.  ENERGIES
C     ARE IN KEV, JQ = SPECIES INDEX = IONIC CHARGE + 1, AND
C     JN = SHELL INDEX.
C
C     USES THE R. MORE REVISED SHIELDING CONSTANTS (AND NEW
C     CALCULATION PROCEDURE) IN LIEU OF THE ORIGINAL 'XSNQ'
C     TREATMENT BASED ON THE H. MAYER SHIELDING CONSTANTS.
C     (ORIGINAL TREATMENT, REFS. 1-4; PRESENT TREATMENT, REF. 5).
C     CERTAIN ASPECTS OF THE FORMALISM (E.G. QN AND EN ARRAYS) ARE
C     'HISTORICAL' IN NATURE.
C
C
C          SIGMA(JN,JM)   R. MORE SHIELDING CONSTANTS
C
C          QN(JQ,JN)      EFFECTIVE NUCLEAR CHARGE AT SHELL JN
C                         CONSIDERING ALL ELECTRONS IN THE GROUND STATE
C                         ION, NOW COMPUTED FROM ZSEN (RATHER THAN
C                         THAN VICE-VERSA).  THIS IS *NOT* THE
C                         SAME AS QN IN MORE'S NOTATION, WHICH
C                         IS THE SHIELDED CHARGE DUE TO INNER
C                         SHIELDING ONLY.
C
C          EN(JQ,JN)      ENERGY LEVELS FOR RECOMBINATION OF THE
C                         GROUND STATE ION.
C
C          EIN(JQ,JN)     IONIZATION POTENTIAL FOR ELECTRON IN SHELL
C                         JN CALCULATED FROM THE DIFFERENCE IN TOTAL
C                         ION ENERGIES WITH AND WITHOUT AN ELECTRON
C                         REMOVED FROM SHELL JN.
C
C          ENM(JQ,JN,JM)  EXCITATION ENERGY FROM SHELL JN UP TO SHELL
C                         JM CALULATED FROM THE DIFFERENCE IN
C                         THE TOTAL ION ENERGIES OF OF THE GROUND
C                         STATE AND EXCITED STATE ELECTRON POPULATIONS.
C
C          ZBETOT(JQ)     TOTAL BINDING ENERGY FOR EACH ION.
C
C
C     REFERENCES:
C
C     1) POST ET AL., PPPL-1352 (CORRECTED)
C     2) LOKKE AND GRASBERGER, LLL REPORT UCRL-52276 (CORRECTED)
C     3) H. MAYER, 'METHODS OF OPACITY CALCULATIONS', LOS ALAMOS
C        SCIENTIFIC  LABORATORY LA-647 (1947).
C     4) XSNQ LISTING
C     5) R. MORE 'ATOMIC PHYSICS IN INERTIAL CONFINEMENT FUSION -
C        PART I' LLL REPORT UCRL-84991 (1981) (TO APPEAR
C        IN 'APPLIED ATOMIC COLLISION PHYSICS', VOLUME II,
C        ACADEMIC PRESS)
C
C     NOTE: THERE ARE DISCREPANCIES BETWEEN THESE REFERENCES FOR A FEW
C           OF THE H. MAYER SHIELDING CONSTANT VALUES.
C
C
      DIMENSION ZSEN(10), ZBETOT(100), ZDUM(10)
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
C
C     MAIN SPECIES LOOP
C
      DO 100 JQ = 1, NSPC
      JJQ = JQ
      IVALNC = NVALNC(JQ)
C
C     COMPUTE ZBETOT AND ZSEN
C
      IHOLE = 0
      IEXCIT = 0
      CALL ADSHEN(JJQ, IHOLE, IEXCIT, ZSEN(1), ZBETOT(JJQ))
C
C     COMPUTE EIN, ENM, AND EN FROM DIFFERENCES IN TOTAL
C     ION ENERGIES (RATHER FROM ZSEN'S) IN ORDER TO AGREE
C     WITH RESULTS IN MORE'S FIGURE III-7.  COMPUTE
C     QN FROM THE ZSEN IN ORDER TO HAVE THE APPROPRIATE LIMITING
C     INTEGRAL CHARGES AT HIGH N SHELLS.
C
C     CLEAR ARRAYS
C
      DO 30 JN = 1, 5
C
      EIN(JQ,JN) = 0.0
C
      DO 20 JM = 1, 10
      ENM(JQ,JN,JM) = 0.0
   20 CONTINUE
   30 CONTINUE
C
C     SKIP FULLY STRIPPED STATE
C
      IF(JQ. EQ. NSPC) GO TO 60
C
      DO 50 JN = 1, IVALNC
C
      IHOLE = JN
      IEXCIT = 0
      CALL ADSHEN(JJQ, IHOLE, IEXCIT, ZDUM(1), ZBEIZ)
      EIN(JQ,JN) = ZBETOT(JQ) - ZBEIZ
C
      IIJM = IVALNC
      IF(JN .EQ. IVALNC .OR. AQN(JQ,IVALNC) .EQ. 0.0) IIJM = IIJM+1
C
      DO 40 JM = IIJM, 10
      IHOLE = JN
      IEXCIT = JM
      CALL ADSHEN(JJQ, IHOLE, IEXCIT, ZDUM(1), ZBEEX)
      ENM(JQ,JN,JM) = ZBETOT(JQ) - ZBEEX
   40 CONTINUE
   50 CONTINUE
C
C     CLEAR ARRAYS
C
   60 DO 70 JN = 1, 10
      EN(JQ,JN) = 0.0
      QN(JQ,JN) = 0.0
   70 CONTINUE
C
C     SKIP NEUTRAL STATE
C
      IF(JQ .EQ. 1) GO TO 100
C
C     COMPUTE EN AND QN
C
      IIJN = IVALNC
      IF(AQN(JQ,IVALNC) .EQ. 0.0) IIJN = IVALNC+1
C
      DO 80 JN = IIJN, 10
      IHOLE = 0
      IEXCIT = JN
      CALL ADSHEN(JJQ, IHOLE, IEXCIT, ZDUM(1), ZBEREC)
      EN(JQ,JN) = ZBEREC - ZBETOT(JQ)
      QN(JQ,JN) = FLOAT(JN) * SQRT(ZSEN(JN)/0.0136)
   80 CONTINUE

  100 CONTINUE
C
      RETURN
      END
C
C---------------------------------------------------------------
C
      SUBROUTINE ADSHEN(KQ, KHOLE, KEXCIT, PSEN, PBETOT)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 1 SEPTEMBER 1981
C
C     INPUT: NUCZ, KQ, KHOLE, KEXCIT, APN
C     OUTPUT: PSEN(10), PBETOT
C
C     COMPUTE SHELL ENERGIES USING NEW SHIELDING CONSTANTS
C     AND METHOD DUE TO R. MORE (UCRL-84991).
C
C     NUCZ = ION NUCLEAR CHARGE
C     KQ = SPECIES INDEX = IONIC CHARGE + 1
C     KHOLE, KEXCIT = ELECTRON REMOVED FROM SHELL KHOLE AND/OR
C                     ADDED (EXCITED) TO SHELL KEXCIT BEFORE
C                     ENERGIES ARE COMPUTED (KHOLE = 0 NO ELECTRON
C                     REMOVED, KEXCIT = 0 NO ELECTRON ADDED)
C     APN(100,10) = ARRAY OF SHELL POPULATIONS
C     PSEN = OUTPUT ARRAY OF SHELL ENERGIES (KEV) (DEFINED POSITIVE)
C     PBETOT = OUTPUT TOTAL ION ENERGY (DEFINED POSITIVE)
C
      DIMENSION PSEN(10), ZAPN(10), ZQ(10), ZEN0(10)
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / SMORE / SIGMOR(10,10)
C
      DATA ((SIGMOR(JN,JM),JN=1,10),JM=1,10)/
     1 0.3125, 0.9380, 0.9840, 0.9954, 0.9970,
     1 0.9970, 0.9990, 0.9999, 0.9999, 0.9999,
     2 0.2345, 0.6038, 0.9040, 0.9722, 0.9979,
     2 0.9880, 0.9900, 0.9990, 0.9990, 0.9999,
     3 0.1093, 0.4018, 0.6800, 0.9155, 0.9796,
     3 0.9820, 0.9860, 0.9900, 0.9920, 0.9999,
     4 0.0622, 0.2430, 0.5150, 0.7100, 0.9200,
     4 0.9600, 0.9750, 0.9830, 0.9860, 0.9900,
     5 0.0399, 0.1597, 0.3527, 0.5888, 0.7320,
     5 0.8300, 0.9000, 0.9500, 0.9700, 0.9800,
     6 0.0277, 0.1098, 0.2455, 0.4267, 0.5764,
     6 0.7248, 0.8300, 0.9000, 0.9500, 0.9700,
     7 0.0204, 0.0808, 0.1811, 0.3184, 0.4592,
     7 0.6098, 0.7374, 0.8300, 0.9000, 0.9500,
     8 0.0156, 0.0624, 0.1392, 0.2457, 0.3711,
     8 0.5062, 0.6355, 0.7441, 0.8300, 0.9000,
     9 0.0123, 0.0493, 0.1102, 0.1948, 0.2994,
     9 0.4222, 0.5444, 0.6558, 0.7553, 0.8300, 0.0100, 0.0400,
     & 0.0900, 0.1584, 0.2450, 0.3492, 0.4655, 0.5760, 0.6723, 0.7612/
C
      DO 10 JN = 1, 10
      ZAPN(JN) = APN(KQ,JN)
   10 CONTINUE
C
C     REMOVE ELECTRON FROM SHELL KHOLE
C
      IF(KHOLE .GT. 0 .AND. ZAPN(KHOLE) .GT. 0.0)
     &         ZAPN(KHOLE) = ZAPN(KHOLE) - 1.0
C
C     ADD ELECTRON TO SHELL KEXCIT
C
      IF(KEXCIT .GT. 0) ZAPN(KEXCIT) = ZAPN(KEXCIT) + 1.0
C
C     COMPUTE NUCLEAR CHARGE AS SHIELDED BY INNER ELECTRONS ONLY
C
      DO 30 JN = 1, 10
C
      ZQ(JN) = FLOAT(NUCZ) - 0.5 * SIGMOR(JN,JN) * ZAPN(JN)
C
      IF(JN .EQ. 1) GO TO 30
      JNM1 = JN - 1
C
      DO 20 JM = 1, JNM1
      ZQ(JN) = ZQ(JN) - SIGMOR(JN,JM) * ZAPN(JM)
   20 CONTINUE
   30 CONTINUE
C
C     COMPUTE POTENTIAL DUE TO OUTER ELECTRONS
C
      DO 50 JN = 1, 10
      ZEN0(JN) = 0.0136 * SIGMOR(JN,JN) * ZAPN(JN)
     &           * ZQ(JN) / (FLOAT(JN)**2)
C
      IF(JN .EQ. 10) GO TO 50
      JNP1 = JN + 1
C
      DO 40 JM = JNP1, 10
      ZEN0(JN) = ZEN0(JN) + 0.0272 * ZAPN(JM) * SIGMOR(JM,JN)
     &                      * ZQ(JM) / (FLOAT(JM)**2)
   40 CONTINUE
C
   50 CONTINUE
C
C     COMPUTE FINAL SHELL ENERGIES AND TOTAL ION BINDING ENERGY,
C     BOTH DEFINED POSITIVE.
C
      PBETOT = 0.0
C
      DO 60 JN = 1, 10
      PSEN(JN) = 0.0136 * ((ZQ(JN)/FLOAT(JN))**2) - ZEN0(JN)
      PBETOT = PBETOT + 0.0136 * ZAPN(JN) * ((ZQ(JN)/FLOAT(JN))**2)
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADMAYR
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     INPUT : NUCZ, NSPC, APN
C     OUTPUT : SIGMA, QN, EN, EIN, ENM
C
C     COMPUTES ENERGY LEVELS FOR THE GROUND STATE IONIC SPECIES
C     AND FOR SPECIES WITH SINGLE INNER ELECTRON VACANCIES.
C     EXCITATION ENERGIES AND IONIZATION POTENTIALS ARE THUS
C     TABULATED FOR ALL OCCUPIED SHELLS FOR EACH SPECIES.  ENERGIES
C     ARE IN KEV, JQ = SPECIES INDEX = IONIC CHARGE + 1, AND
C     JN = SHELL INDEX.
C
C     FOR EACH CONFIGURATION AND SHELL OF INTEREST, THE EFFECTIVE
C     NUCLEAR CHARGE IS CALCULATED USING THE H. MAYER SHIELDING
C     CONSTANTS.  THE ENERGY LEVELS ARE THEN FOUND USING THIS EFFECTIVE
C     CHARGE IN THE PAULI APPROXIMATION TO THE DIRAC ENERGY LEVEL
C     FORMULA.  THE SHELL ENERGY IS TAKEN AS THE MEAN OF THE LOWEST
C     AND HIGHEST SUBSHELL ENERGIES FOR THAT N.
C
C
C          SIGMA(JN,JM)   H. MAYER SHIELDING CONSTANTS
C
C          QN(JQ,JN)      EFFECTIVE NUCLEAR CHARGE AT SHELL JN
C                         CONSIDERING ALL ELECTRONS IN THE GROUND STATE ION.
C
C          EN(JQ,JN)      ENERGY LEVELS FOR GROUND STATE ION
C                         (FORMERLY EQN(JQ,JN) IN PREVIOUS VERSIONS)
C
C          EIN(JQ,JN)     IONIZATION POTENTIAL FOR ELECTRON IN SHELL JN;
C                         USES SHIELDED CHARGE CALCULATED WITH ONE
C                         ELECTRON REMOVED FROM SHELL JN.
C
C          ENM(JQ,JN,JM)  EXCITATION ENERGY FROM SHELL JN UP TO SHELL JM;
C                         THE DIFFERENCE BETWEEN THE JN AND JM ENERGY LEVELS
C                         CALCULATED WITH THE JN SHELL ELECTRON REMOVED.
C
C
C
C     REFERENCES:
C
C     1) POST ET AL., PPPL-1352 (CORRECTED)
C     2) LOKKE AND GRASBERGER, LLL REPORT UCRL-52276 (CORRECTED)
C     3) H. MAYER, 'METHODS OF OPACITY CALCULATIONS', LOS ALAMOS
C        SCIENTIFIC  LABORATORY LA-647 (1947).
C     4) XSNQ LISTING
C
C     NOTE: THERE ARE DISCREPANCIES BETWEEN THESE REFERENCES FOR A FEW
C           SHIELDING CONSTANT VALUES.
C
C
      DIMENSION ZSIGMA(10,10)
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
C
      COMMON / SMAYER / SIGMAY(10,10)
C
      DATA ((SIGMAY(JN,JM),JN=1,10),JM=1,10)/
     1 .625,.938,.981,.987,.994,.997,.999,1.00,1.00,1.00,
     2 .235,.690,.893,.940,.970,.984,.990,.993,.995,1.00,
     3 .109,.397,.702,.850,.920,.955,.970,.980,.990,1.00,
     4 .0617,.235,.478,.705,.830,.900,.950,.970,.980,.990,
     5 .0398,.155,.331,.531,.720,.830,.900,.950,.970,.980,
     6 .0277,.109,.239,.400,.580,.735,.830,.900,.950,.970,
     7 .0204,.0808,.178,.310,.459,.610,.745,.830,.900,.950,
     8 .0156,.0625,.138,.243,.371,.506,.635,.750,.830,.900,
     9 .0123,.0494,.111,.194,.299,.431,.544,.656,.760,.830,
     & .010,.040,.090,.158,.245,.353,.460,.576,.670,.765/
C
C
C
C     START MAIN SPECIES LOOP
C
C
      DO 100 JQ = 1 , NSPC
C
      IVALNC = NVALNC(JQ)
C
C
C     COMPUTE QN AND EN
C
      DO 30 JN = 1 , 10
      ZJN = FLOAT(JN)
C
      ZQN = FLOAT(NUCZ)
C
      DO 20 JM = 1 , IVALNC
      ZQN = ZQN - SIGMAY(JN,JM) * APN(JQ,JM)
   20 CONTINUE
C
      ZQN = ZQN + APN(JQ,JN) * SIGMAY(JN,JN) / (2. * ZJN**2)
C
      QN(JQ,JN) = ZQN
      IN = JN
      EN(JQ,JN) = ADPALI(ZQN,IN)
   30 CONTINUE
C
C
C     CLEAR ARRAYS THEN COMPUTE EIN AND ENM
C
C
      DO 40 JN = 1 , 5
      EIN(JQ,JN) = 0.
      DO 40 JM = 1 , 10
      ENM(JQ,JN,JM) = 0.
   40 CONTINUE
C
C     SKIP IONIZATION AND EXCITATION ENERGIES FOR FULLY STRIPPED ION
C
      IF(JQ .EQ. NSPC) GO TO 100
C
C
      DO 70 JN = 1 , IVALNC
      ZJN = FLOAT(JN)
C
C     COMPUTE ZQNN = EFFECTIVE CHARGE AT SHELL JN WITH ONE ELECTRON
C     REMOVED FROM SHELL JN.
C
      ZQNN = QN(JQ,JN) + SIGMAY(JN,JN) * (1. - 1./(2.*ZJN**2))
C
C     EIN(JQ,JN) = IONIZATION POTENTIAL FOR SHELL JN
C
      IN = JN
      EIN(JQ,JN) = ADPALI(ZQNN,IN)
C
C     COMPUTE ENM(JQ,JN,JM) FROM ZQNM, EFFECTIVE CHARGE AT SHELL
C     JM WITH ONE ELECTRON REMOVED FROM SHELL JN.
C
      IIJM = IVALNC
      IF(JN .EQ. IVALNC .OR.
     &         AQN(JQ,IVALNC) .EQ. 0.) IIJM = IVALNC + 1
C
      DO 60 JM = IIJM , 10
C
      ZQNM = QN(JQ,JM) + SIGMAY(JM,JN)
C
      IM = JM
      ENM(JQ,JN,JM) = EIN(JQ,JN) - ADPALI(ZQNM,IM)
C
   60 CONTINUE
C
C
   70 CONTINUE
C
C
  100 CONTINUE
C
C
      RETURN
      END
C
C
C------------------------------------------------------------------
C
C
      FUNCTION ADPALI(PQ,KN)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     EVALUATE THE PAULI APPROXIMATION TO THE
C     THE DIRAC ENERGY LEVEL EQUATION, TAKING THE MEAN OF THE
C     LOWEST AND HIGHEST SUBSHELL ENERGIES TO BE THE SHELL ENERGY.
C     ADPALI(PQ,KN) IS IN KEV,  PQ IS THE EFFECTIVE NUCLEAR
C     CHARGE (IN UNITS OF ELECTRONIC CHARGE) AND KN IS THE PRINCIPAL
C     QUANTUM NUMBER.
C
C
      ZALPHA = .007298
C
C
      ZKN = FLOAT(KN)
C
      ADPALI = .0136 * ( (PQ / ZKN)**2 ) *
     &     (1. + ((ZALPHA*PQ)**2) * (0.5 - 1./(4.*ZKN)) / ZKN)
C
C
      RETURN
      END
      SUBROUTINE ADTIP(KTIPOK)
C
C     SUBSTITUTE TABULATED IONIZATION POTENTIALS FOR THOSE
C     CALCULATED BY SHIELDING CONSTANTS IN SETENM.  NOTE THAT THE
C     'NUMBER OF EQUIVALENT ELECTRONS' IS NOT CORRECTLY RE-DEFINED
C     TO ACCOMPANY THIS PROCEDURE, BUT THE APPROXIMATION IS STILL
C     AN IMPROVEMENT FOR NEAR-NEUTRAL SPECIES.
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      DIMENSION ZTIP(100)
C
      DO 1 JSPC = 1, 100
      ZTIP(JSPC) = 0.0
    1 CONTINUE
C
      IF(NUCZ .NE. 2) GO TO 3
C
C     HELIUM
C
      ZTIP(1) = 24.58E-3
      ZTIP(2) = 54.4E-3
C
      GO TO 100
    3 IF(NUCZ .NE. 3) GO TO 4
C
C     LITHIUM
C
      ZTIP(1) = 5.4E-3
      ZTIP(2) = 75.6E-3
      ZTIP(3) = 122.4E-3
C
      GO TO 100
C
    4 IF(NUCZ .NE. 4) GO TO 5
C
C     BERYLLIUM
C
      ZTIP(1) = 9.32E-3
      ZTIP(2) = 18.2E-3
      ZTIP(3) = 153.9E-3
      ZTIP(4) = 217.6E-3
C
      GO TO 100
C
    5 IF(NUCZ .NE. 5) GO TO 6
C
C     BORON
C
      ZTIP(1) = 8.296E-3
      ZTIP(2) = 25.15E-3
      ZTIP(3) = 37.9E-3
      ZTIP(4) = 259.3E-3
      ZTIP(5) = 340.0E-3
C
      GO TO 100
C
    6 IF(NUCZ .NE. 6) GO TO 8
C
C     CARBON
C
      ZTIP(1) = 11.3E-3
      ZTIP(2) = 24.4E-3
      ZTIP(3) = 47.9E-3
      ZTIP(4) = 64.5E-3
      ZTIP(5) = 0.392
      ZTIP(6) = 0.490
C
      GO TO 100
    8 IF(NUCZ .NE. 8) GO TO 10
C
C     OXYGEN
C
      ZTIP(1) = 13.6E-3
      ZTIP(2) = 35.1E-3
      ZTIP(3) = 54.9E-3
      ZTIP(4) = 77.4E-3
      ZTIP(5) = 0.114
      ZTIP(6) = 0.138
      ZTIP(7) = 0.739
      ZTIP(8) = 0.871
C
      GO TO 100
   10 IF(NUCZ .NE. 10) GO TO 13
C
C     NEON
C
      ZTIP(1) = 21.6E-3
      ZTIP(2) = 41.0E-3
      ZTIP(3) = 63.5E-3
      ZTIP(4) = 97.1E-3
      ZTIP(5) = 0.126
      ZTIP(6) = 0.158
      ZTIP(7) = 0.207
      ZTIP(8) = 0.239
      ZTIP(9) = 1.196
      ZTIP(10) = 1.362
C
      GO TO 100
   13 IF(NUCZ .NE. 13) GO TO 18
C
C     ALUMINUM
C
      ZTIP(1) = 5.99E-3
      ZTIP(2) = 18.8E-3
      ZTIP(3) = 28.4E-3
      ZTIP(4) = 0.120
      ZTIP(5) = 0.154
      ZTIP(6) = 0.190
      ZTIP(7) = 0.241
      ZTIP(8) = 0.285
      ZTIP(9) = 0.330
      ZTIP(10) = 0.399
      ZTIP(11) = 0.442
      ZTIP(12) = 2.086
      ZTIP(13) = 2.304
C
      GO TO 100
   18 IF(NUCZ .NE. 18) GO TO 22
C
C     ARGON
C
      ZTIP(1) = 15.8E-3
      ZTIP(2) = 27.6E-3
      ZTIP(3) = 40.7E-3
      ZTIP(4) = 59.8E-3
      ZTIP(5) = 75.0E-3
      ZTIP(6) = 91.0E-3
      ZTIP(7) = 0.124
      ZTIP(8) = 0.143
      ZTIP(9) = 0.422
      ZTIP(10) = 0.479
      ZTIP(11) = 0.539
      ZTIP(12) = 0.618
      ZTIP(13) = 0.686
      ZTIP(14) = 0.756
      ZTIP(15) = 0.855
      ZTIP(16) = 0.918
      ZTIP(17) = 4.121
      ZTIP(18) = 4.426
C
      GO TO 100
   22 IF(NUCZ .NE. 22) GO TO 26
C
C     TITANIUM
C
      ZTIP(1) = 6.82E-3
      ZTIP(2) = 13.6E-3
      ZTIP(3) = 27.5E-3
      ZTIP(4) = 43.3E-3
      ZTIP(5) = 99.2E-3
      ZTIP(6) = 0.119
      ZTIP(7) = 0.141
      ZTIP(8) = 0.169
      ZTIP(9) = 0.193
      ZTIP(10) = 0.216
      ZTIP(11) = 0.265
      ZTIP(12) = 0.292
      ZTIP(13) = 0.787
      ZTIP(14) = 0.861
      ZTIP(15) = 0.940
      ZTIP(16) = 1.042
      ZTIP(17) = 1.131
      ZTIP(18) = 1.220
      ZTIP(19) = 1.342
      ZTIP(20) = 1.425
      ZTIP(21) = 6.249
      ZTIP(22) = 6.626
C
      GO TO 100
   26 IF(NUCZ .NE. 26) GO TO 100
C
C     IRON
C
      ZTIP(1) = 7.87E-3
      ZTIP(2) = 16.2E-3
      ZTIP(3) = 30.7E-3
      ZTIP(4) = 54.8E-3
      ZTIP(5) = 75.0E-3
      ZTIP(6) = 99.0E-3
      ZTIP(7) = 0.125
      ZTIP(8) = 0.151
      ZTIP(9) = 0.235
      ZTIP(10) = 0.262
      ZTIP(11) = 0.290
      ZTIP(12) = 0.331
      ZTIP(13) = 0.361
      ZTIP(14) = 0.392
      ZTIP(15) = 0.457
      ZTIP(16) = 0.490
      ZTIP(17) = 1.265
      ZTIP(18) = 1.358
      ZTIP(19) = 1.456
      ZTIP(20) = 1.582
      ZTIP(21) = 1.689
      ZTIP(22) = 1.799
      ZTIP(23) = 1.950
      ZTIP(24) = 2.045
      ZTIP(25) = 8.828
      ZTIP(26) = 9.278
C
C     END OF CURRENT TABULATION; RETURN IF SELECTED
C     ELEMENT NOT PRESENT AND SET FLAG.
C
  100 KTIPOK = 0
      IF(ZTIP(1) .EQ. 0.0) RETURN
C
      KTIPOK = 1
C
      DO 150 JSPC = 1, NUCZ
      IVALNC = NVALNC(JSPC)
      JSPCP1 = JSPC + 1
C
      EIN(JSPC,IVALNC) = ZTIP(JSPC)
      EN(JSPCP1,IVALNC) = ZTIP(JSPC)
  150 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADFNM
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 8/10/79  CORRECT ZEPSNN(59) (4F13) VALUE FROM 3.43E-3
C                   TO 3.43E-2 (ERROR IN POST ET AL.)
C     REV: 6/15/78
C
C     INPUT : NUCZ, APN, AQN, NVALNC, ENM
C     OUTPUT : FNM, FNN, ENN
C
C     CALCULATES ABSORPTION OSCILLATOR STRENGTHS FOR TRANSITIONS
C     FROM SHELL IN TO SHELL IM FOR EACH SPECIES JQ AND STORES THEM
C     IN THE ARRAY FNM(JQ,IN,IM).  F(N,M) = (HYDROGENIC VALUE,
C     TABULATED FOR LOW N,M CALCULATED OTHERWISE) X (CORRECTION
C     FACTOR(S)) X (ELECTRON POPULATION OF SHELL IN) X (VACANY
C     FRACTION OF SHELL IM).
C
C     FOR N - N (E.G. DELTA N EQUAL ZERO) TRANSITIONS WE HAVE
C     F(N,N) = ANN + BNN/ZNUC, WHICH REPRESENTS A FIT TO THE
C     'EQUIVALENT' SINGLE TRANSITION USED TO MODEL ALL THE PHYSICAL
C     N - N TRANSITIONS.  ALSO COMPUTED IS E(N,N), THE 'EQUIVALENT'
C     ENERGY, E(N,N) = EPSILONNN * (ZCHG+1) ** ALPHANN.  THE ZANN, ETC.
C     ARRAYS ARE INDEXED BY THE TOTAL NUMBER OF BOUND ELECTRONS
C     IN THEIR ELECTRONIC CONFIGURATIONS (SEE REFERENCE (1)).
C
C     NOTE THAT HERE, AS IN OTHER ROUTINES, THE FILLING OF SHELLS
C     IN SIMPLE ASCENDING ORDER IS ASSUMED BOTH EXPLICITLY
C     (E.G., IN THE APN ARRAY CONTENTS) AS WELL AS IMPLICITY (THROUGH
C     THE USE OF NVALNC, ETC.).  ALSO, THE ANN, BNN, ETC.
C     COEFFICIENTS ARE DEFINED SPECIFICALLY FOR SUCH STATES.
C
C
C     REFERENCES:
C     1) POST ET AL., PPPL-1352 (WITH CORRECTIONS)
C     2) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #4, 7/29/76
C     3) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #10, 1/5/77
C     4) LOKKE & GRASBERGER, 'XSNQ-U A NON-LTE EMISSION AND
C        ABSORPTION COEFFICIENT SUBROUTINE', LLL REPORT UCRL-52276
C     5) MERTS ET AL., 'THE CALCULATED POWER OUTPUT FROM
C        A THIN IRON-SEEDED PLASMA', LOS ALAMOS REPORT LA-6220-MS
C     6) XSNQ FORTRAN LISTING
C
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
C
      COMMON / FNMDAT / ZANN(60), ZBNN(60), ZEPSNN(60), ZALPNN(60),
     &          ZY0(10), ZY1(10), ZY2(10)
C
C
      DATA (ZANN(J),J=1,60) / 0., 0.,
     1 -.02, -.01, -.04, .04, -.01, -.01, .01, 0.,
     2  .01, .1, .25, .17, .08, -.02, -.14, -.27,
     3 -.29, -.3, -.3, -.29, -.27, -.24, -.2, -.14,
     4 -.08, 0., .97, 1.96, 1.92, 1.89, 1.86, 1.83,
     5 1.78, 1.73, 1.41, 1.05, .674, .261, -.167,
     6 -.637, -1.14, -1.67, -2.26, -2.88, -2.9,
     7 -2.83, -2.72, -2.61, -2.45, -2.27, -2.05, -1.81,
     8 -1.55, -1.25, -.94, -.63, -.31, 0. /
C
C
      DATA (ZBNN(J),J=1,60) / 0., 0.,
     1 2., 4.33, 4.68, 3.6, 2.85, 2.4, 1.3, 0.,
     2 10.5, 20., 16.4, 24.7, 34.1, 44.8, 56.8, 70.3,
     3 68.6, 65.8, 61.9, 56.9, 50.7, 45.5, 34.4,
     4 24.2, 12.8, 0., -25.8, -53.1, -28.1, -2.97,
     5 23., 50.3, 79.1, 109.2, 141.7, 176., 213.9,
     6 256.6, 299.6, 348.2, 400.5, 456.6, 518.1,
     7 584.0, 571.9, 550.7, 524.9, 496.8, 463.4,
     8 427., 385.3, 339.8, 291.3, 237.4, 181.3, 122.9,
     9 62.2, 0. /
C
C
      DATA (ZEPSNN(J),J=1,28) / 0., 0.,
     1 2.04E-3, 4.49E-3, 6.8E-3, 8.16E-3, 6.80E-3,
     2 1.06E-2, 1.31E-2, 0., 2.3E-3, 3.88E-3,
     3 5.71E-3, 5.44E-3, 8.16E-3, 6.8E-3, 8.3E-3,
     4 4.32E-3, 5.11E-3, 6.04E-3, 7.08E-3,
     5 8.12E-3, 9.69E-3, 1.13E-2, 1.37E-2, 1.49E-2,
     6 1.71E-2, 0. /
C
      DATA (ZEPSNN(J),J=29,60) / 1.52E-5, 1.93E-5,
     1 5.62E-5, 1.20E-4, 2.13E-4, 3.34E-4, 4.66E-4,
     2 6.66E-4, 9.10E-4, 1.25E-3, 1.69E-3, 1.98E-3,
     3 3.03E-3, 3.92E-3, 5.19E-3, 6.40E-3, 8.05E-3,
     4 9.80E-3, 1.08E-2, 1.19E-2, 1.32E-2, 1.47E-2,
     5 1.59E-2, 1.81E-2, 1.92E-2, 2.10E-2, 2.33E-2,
     6 2.57E-2, 2.82E-2, 3.10E-2, 3.43E-2, 0. /
C
C
      DATA (ZALPNN(J),J=1,60) / 1.0, 1.0,
     1 1.0, .93, .86, .80, .87, .77, .77, 1.0,
     2 .99, .90, .83, .87, .77, .82, .79, 1.06,
     3 1.03, 1., .97, .94, .91, .88, .85,
     4 .83, .80, 1.0, 2.46, 2.4, 2.14, 1.96,
     5 1.83, 1.72, 1.64, 1.57, 1.49, 1.41, 1.34,
     6 1.30, 1.19, 1.13, 1.06, 1.01, .95, .91,
     7 .89, .87, .85, .83, .82, .80, .78, .76,
     8 .74, .72, .70, .68, .66, 1.0 /
C
C
      DATA (ZY0(J),J=1,10) / 1.0, .511, .400, 7*.230 /
C
      DATA (ZY1(J),J=1,10) / 0., .116, .0228, 7*.00844 /
C
      DATA (ZY2(J),J=1,10) / 0., -.00689, .000596, 7*.000488 /
C
C
C     *********************************************************
C
C
C     ZERO OSCILLATOR STRENGTH ARRAYS AND ENN ARRAY
C
      DO 1 JQ = 1 , 100
      FNN(JQ) = 0.
      ENN(JQ) = 0.
      DO 1 JN = 1 , 5
      DO 1 JM = 1 , 10
    1 FNM(JQ,JN,JM) = 0.
C
      ZNUC = FLOAT ( NUCZ )
C
C
C     ******************************************************
C     ***** LOOP FOR ALL SPECIES EXCEPT FULLY STRIPPED *****
C     ******************************************************
C
C
      DO 100 JQ = 1 , NUCZ
C
      ZCHG = FLOAT(JQ-1)
      IBOUND = NUCZ - JQ + 1
      IVALNC = NVALNC(JQ)
C
C
C     ********** CALCULATE F(N,M) ( M > N ) **********
C
C
C     BEGIN BY EVALUTING HYDROGENIC F(N,M) EQUATION, WITH
C     EXPLICIT VALUES SUBSTITUTED FOR N < 5, M <= 5.
C
C
      DO 10 JN = 1 , IVALNC
      IIJM = MAX0 ( JN + 1 , IVALNC )
      ZN = FLOAT(JN)
C
      DO 10 JM = IIJM , 10
      ZM = FLOAT(JM)
      FNM(JQ,JN,JM) = 1.96 * ZN * ( ( ZM/(ZM*ZM - ZN*ZN) )**3 )
   10 CONTINUE
C
C     SUBSTITUTE TABLED VALUES FOR LOW N , M AS PER REF #4
C
      FNM(JQ,1,2) = 0.4161
      FNM(JQ,1,3) = 0.0792
      FNM(JQ,1,4) = 0.0290
      FNM(JQ,1,5) = 0.0139
      FNM(JQ,2,3) = 0.6370
      FNM(JQ,2,4) = 0.1190
      FNM(JQ,2,5) = 0.0443
      FNM(JQ,3,4) = 0.8408
      FNM(JQ,3,5) = 0.1499
      FNM(JQ,4,5) = 1.037
C
C     MODIFY NONZERO F(N,N+1) VALUES FOR N > 1 AS PER REF #3
C
      DO 20 J = 1 , 2
      JN = IVALNC + J - 2
      IF(JN .LE. 1) GO TO 20
      JM = JN + 1
      IF(AQN(JQ,JM) .EQ. 0.) GO TO 20
C
      ZNCHG = ZCHG + APN(JQ,JM)
      ZEHYDR = .0136 * ((ZNCHG+1.)**2) *
     &                   ( 1./FLOAT(JN)**2 - 1./FLOAT(JM)**2 )
      ZRATIO = ZEHYDR / ENM(JQ,JN,JM)
      IF(ZRATIO .GT. 1.0) ZRATIO = 1.0
      ZYNPN = ZY0(JN) + ZY1(JN) * APN(JQ,JN) + ZY2(JN) * APN(JQ,JN)**2
C
      FNM(JQ,JN,JM) = FNM(JQ,JN,JM) * ZRATIO * ZYNPN
C
   20 CONTINUE
C
C
C     MULTIPLY BY ELECTRON POPULATION OF SHELL N AND FRACTIONAL
C     VACANCY OF SHELL M
C
C
      DO 30 JN = 1 , 5
      DO 30 JM = 1 , 10
      FNM(JQ,JN,JM) = APN(JQ,JN) * AQN(JQ,JM) * FNM(JQ,JN,JM)
   30 CONTINUE
C
C
C
C     ***** CALCULATE N - N 'EFFECTIVE' OSCILLATOR STRENGTHS *****
C                   ***** AND TRANSITION ENERGIES *****
C
C
C
      INDEX = IBOUND
C     IF MORE THAN 60 BOUND ELECTRONS USE CONSTANTS FOR ION WITH
C     32 FEWER ELECTRONS
      IF(INDEX.GT.60)INDEX = INDEX - 32
      IF(INDEX.GT.60)INDEX = INDEX - 32
C
      FNN(JQ) = ZANN(INDEX) + ZBNN(INDEX) / ZNUC
C
      ENN(JQ) = ZEPSNN(INDEX) * (ZCHG + 1.)**ZALPNN(INDEX)
C
C
C
C
C     CHECK TO INSURE ALL F VALUES ARE POSITIVE
C
      IF(FNN(JQ) .GE. 0.) GO TO 40
      FNN(JQ) = 0.
      WRITE(6,900) JQ , IVALNC , IVALNC
  900 FORMAT(1X,26H***** NEGATIVE F VALUE AT ,3I4,6H *****)
   40 DO 41 JN = 1 , 5
      DO 41 JM = 1 , 10
      IF(FNM(JQ,JN,JM) .GE. 0.) GO TO 41
      FNM(JQ,JN,JM) = 0.
      WRITE(6,900) JQ , JN , JM
   41 CONTINUE
C
C     DIAGNOSTIC WRITE STATMENTS TO DISPLAY F VALUES
C
CW    WRITE(6,901)JQ,FNN(JQ),ENN(JQ)
CW  901 FORMAT(/1X,'SETFNM: JQ =',I3,'  FNN =',1PE9.2,'  ENN =',E9.2)
CW    WRITE(6,902)( JN,(FNM(JQ,JN,JM),JM=1,10),JN=1,5)
CW  902 FORMAT(1X,'JN =',I2,10F6.1)
C
C
C     *************** END OF JQ SPECIES LOOP ***************
C
C
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADERC(PTE,PNE)
C
C     CALL ELECTRON RATE COEFFICIENT ROUTINES
C
      CALL ADECI(PTE,PNE)
      CALL ADDREC(PTE,PNE)
      CALL ADRREC(PTE,PNE)
      CALL ADECEX(PTE,PNE)
      CALL ADBREM(PTE,PNE)
C
      RETURN
      END
      SUBROUTINE ADECI(PTE, PNE)
C
C     CALCULATE THE ELECTRON COLLISIONAL IONIZATION RATE USING
C     THE FORMULAS SELECTED BY THE LECI SWITCH.  LECIOK RETURNS
C     THE RATE FORMULAS ACTUALLY USED, WHICH MAY DIFFER FROM
C     LECI IF LECI WAS SET INAPPROPRIATELY FOR THE GIVEN ELEMENT.
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
      COMMON / ADPARM / LADEN, LADTIP, LTIPOK, LECI, LECIOK
C
C     ILMAX = CURRENT MAXIMUM NUMBER OF RATE OPTIONS
C
      ILMAX = 3
C
C     CHECK FOR VALID LECI; RETURN LECI = 1 OPTION IF NOT VALID
C
      IF(LECI .LT. 1 .OR. LECI .GT. ILMAX) GO TO 10
C
      GO TO (10, 20, 30) LECI
C
C     LECI = 1; USE XSNQ DERIVED IONIZATION RATE CALCULATION
C
   10 CALL ADXECI(PTE, PNE)
      LECIOK = 1
      RETURN
C
C     LECI = 2; USE BELFAST GROUP LIGHT ELEMENT RATE FITS
C
   20 CALL ADBFI(NUCZ, PTE, PNE, RCLION, CIZLOS, IFLAG)
      IF(IFLAG .NE. 0) GO TO 10
      LECIOK = 2
      RETURN
C
C     LECI = 3; USE YOUNGER'S RATES FOR AVAILABLE ELEMENTS AND SHELLS,
C               XSNQ RATES OTHERWISE ON A CHARGE STATE-BY-CHARGE STATE
C               BASIS
C
   30 CALL ADXECI(PTE, PNE)
      CALL ADYGR(NUCZ, PTE, PNE, RCLION, CIZLOS, IFLAG)
      LECIOK = 3
      IF(IFLAG .NE. 0) LECIOK = 1
      RETURN
C
C
C
      END
      SUBROUTINE ADXECI(PTE, PNE)
C
C     REV: 3/30/83 ADECI --> ADXECI TO ACCOMODATE ALTERNATIVE RATES
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     COMPUTES THE XSNQ-DERIVED COLLISIONAL IONIZATION RATES RCLION
C     AT THE ELECTRON TEMPERATURE PTE (KEV), WHERE JQ=(IONIC CHARGE+1).
C     THESE RATES HAVE THE UNITS (SEC-1).
C
C     NOTES:
C
C     1) ELECTRONS ARE ASSUMED TO FILL IONIC SHELLS IN SIMPLE ASCENDING
C        ORDER.  THIS IS ASSUMED BOTH EXPLICITLY (THE APN ARRAY)
C        AS WELL AS IMPLICITLY THROUGH THE USE OF NVALNC VALUES
C        IN LOOP INDICES, ETC.
C     2) THROUGH USE OF THE FUNCTION EXPUND, EXP(X) IS SET EQUAL TO
C        ZERO IN THE RATE EQUATIONS FOR VALUES OF X < -45.
C     3) RATES FOR THE SPECIES JQ REFER TO TRANSFORMATIONS OF
C        THAT SPECIES INTO ADJACENT IONIZATION STATES, NOT THE
C        RATES AT WHICH THAT SPECIES IS FORMED.
C
C
C     PRINCIPAL PHYSICS REFERENCES:
C
C     1) POST, JENSEN, TARTER, GRASBERGER, LOKKE, 'STEADY STATE
C        RADIATIVE COOLING RATES FOR LOW-DENSITY HIGH-TEMPERATURE
C        PLASMAS', PPPL-1352 (SUB. TO N.F.) (WITH CORRECTIONS)
C     2) LOKKE AND GRASBERGER, 'XSNQ-U A NON-LTE EMISSION AND
C        ABSORPTION COEFFICIENT SUBROUTINE', LLL REPORT UCRL-52276
C        (WITH CORRECTIONS).
C     3) XSNQ FORTRAN LISTING
C
C
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
C
      ZJLKEV = 1.6021E-16
C
C
      ZTE = PTE
C
      DO 100 JQ = 1 , NUCZ
C
      ZCHG = FLOAT(JQ - 1)
      IVALNC = NVALNC(JQ)
C
C
C
C     ******************** IONIZATION RATE ********************
C
C
C     CALCULATE ELECTRON COLLISIONAL IONIZATION RATES PER SHELL
C     AND SUM FOR ALL POPULATED SHELLS. ALSO COMPUTE
C     IONIZATION ENERGY LOSS RATE CIZLOS.
C
      ZISUM = 0.
      ZESUM = 0.0
C
      DO 20 JN = 1 , IVALNC
      ZJN = JN
      ZXN = EIN(JQ,JN)/ZTE
C
      ZY = ALOG10(ZXN/4.)
      ZFY = .23 + .046 * ZY + .1074 * ZY**2
     &          - .0459 * ZY**3 - .01505 * ZY**4
      ZGAUNT = 12.18 * EXP(-ZJN / (ZJN+5.)) * (1. + .0335*ZJN)
     &         * (1. - .622/(ZCHG+1.) - .0745/(ZCHG+1.)**2)
     &         * ABS(ZFY)
      ZITEMP = APN(JQ,JN) * EXPUND(-ZXN,-45.) *
     &         (1.-EXPUND(-ZXN,-45.)) * ZGAUNT / EIN(JQ,JN)**2
      ZISUM = ZISUM + ZITEMP
      ZESUM = ZESUM + EIN(JQ,JN) * ZITEMP
   20 CONTINUE
C
      ZISUM = ZISUM * SQRT(ZTE)
      ZESUM = ZESUM * SQRT(ZTE)
      IF(ZISUM .LT. 1.E-26) ZISUM = 0.
      IF(ZESUM .LT. 1.E-26) ZESUM = 0.0
      RCLION(JQ) = PNE * 3.44E-11 * ZISUM
      CIZLOS(JQ) = PNE * ZJLKEV * 3.44E-11 * ZESUM
C
C
  100 CONTINUE
C
C
      RCLION(NUCZ + 1) = 0.
      CIZLOS(NUCZ + 1) = 0.0
C
C
      RETURN
      END
      SUBROUTINE ADBFI(KNUCZ, PTE, PNE, PBFI, PBFIL, KFLAG)
C
C     CALCULATE THE ELECTRON COLLISIONAL IONIZATION RATE COEFFICIENTS
C     USING FITTING FORMULAS FROM THE BELFAST GROUP, CULHAM
C     REPORT CLM-R216 (DECEMBER 1981):
C
C          ATOMIC AND MOLECULAR DATA FOR FUSION, PART I
C
C          RECOMMENDED CROSS SECTIONS AND RATES FOR ELECTRON
C          IONISATION OF LIGHT ATOMS AND IONS
C
C          K L BELL, H B GILBODY, J G HUGHES, A E KINGSTON, F J SMITH
C          THE QUEEN'S UNIVERSITY OF BELFAST, BELFAST, N IRELAND
C
C     KNUCZ = ELEMENT ATOMIC NUMBER ( Z=1 TO Z=8 ONLY)
C     PTE = ELECTRON TEMPERATURE (KEV)
C     PNE = ELECTRON DENSITY (CM-3)
C
C     PBFI = IONIZATION RATES FOR THE Z+1 CHARGE STATES (SEC-1)
C     PBFIL = ARRAY OF IONIZATION ENERGY LOSS RATES (WATTS/ION)
C     KFLAG = 0 FOR O.K., = 1 FOR INVALID ELEMENT
C
      COMMON / COMBFI / EH,     AH(6),    AAH,     BTH(3),
     &                  EHE(2), AHE(6,2), AAHE(2), BTHE(3,2),
     &                  ELI(3), ALI(6,3), AALI(3), BTLI(3,3),
     &                  EBE(4), ABE(6,4), AABE(4), BTBE(3,4),
     &                  EB(5),  AB(6,5),  AAB(5),  BTB(3,5),
     &                  EC(6),  AC(6,6),  AAC(6),  BTC(3,6),
     &                  EN(7),  AN(6,7),  AAN(7),  BTN(3,7),
     &                  EO(8),  AO(6,8),  AAO(8),  BTO(3,8)
C
C
      DIMENSION PBFI(1), PBFIL(1)
C
C----------------------------------------------------------------------
C
C     HYDROGEN
C
      DATA EH / 13.60 /
      DATA AH / 2.3742E-8, -3.6866E-9, -1.0366E-8,
     &         -3.8010E-9,  3.4159E-9,  1.6834E-9 /
      DATA AAH / 2.4617E-8 /
      DATA BTH / 9.5986E-8, -9.2463E-7, 3.9973E-6 /
C
C     HELIUM
C
      DATA EHE / 24.6, 54.42 /
      DATA AHE / 1.4999E-8, 5.6656E-10, -6.0821E-9,
     1          -3.5894E-9, 1.5529E-9,   1.3207E-9,
     2           3.4356E-9, -1.6865E-9, -6.9236E-10,
     2           9.7863E-11, 1.5591E-10, 6.2236E-11 /
      DATA AAHE / 3.1373E-8, 3.0772E-9 /
      DATA BTHE / 4.7893E-8, -7.7359E-7, 3.7366E-6,
     2            1.1902E-8, -1.1514E-7, 5.0489E-7 /
C
C     LITHIUM
C
      DATA ELI / 5.39, 75.64, 122.45 /
      DATA ALI / 9.9655E-8, -5.5941E-8, -5.5228E-8,
     1           4.0589E-8, 1.4800E-8, -1.3120E-8,
     2           3.4023E-9, -7.6588E-10, -8.6078E-10,
     2           -8.9748E-10, 4.1661E-10, 3.3188E-10,
     3           1.1786E-9, -8.7637E-10, -9.3373E-11,
     3           2.1137E-10, 1.9017E-11, -4.0679E-11 /
      DATA AALI / 4.5456E-8, 7.3504E-9, 1.9767E-9 /
      DATA BTLI / 2.7800E-7, -1.5830E-6, 5.4650E-6,
     2            5.3459E-10, -5.6387E-8, 2.9577E-7,
     3           -1.0926E-9, 8.8700E-10, 6.0764E-9 /
C
C     BERYLLIUM
C
      DATA EBE / 9.32, 18.21, 153.89, 217.71 /
      DATA ABE / 7.4206E-8, -1.5520E-8, -3.9403E-8,
     1           7.2155E-9, 1.1098E-8, -2.5501E-9,
     2           1.7136E-8, -1.3997E-8, -2.9656E-11,
     2           3.0777E-9, -1.1676E-10, -5.1930E-10,
     3           1.6039E-9, -6.4336E-10, -7.7804E-10,
     3           3.3527E-10, 2.1889E-10, -1.0600E-10,
     4           4.9587E-10, -3.6870E-10, -3.9284E-11,
     4           8.8928E-11, 8.0007E-12, -1.7114E-11 /
      DATA AABE / 2.1732E-7, 2.7494E-8, 2.7873E-9, 8.3163E-10 /
      DATA BTBE / -2.1649E-7, 2.8120E-7, 5.3031E-7,
     2            -1.5900E-8, 2.5911E-8, 2.4061E-8,
     3            -2.4333E-10, -9.7973E-9, 5.2094E-8,
     4            -4.5966E-10, 3.7317E-10, 2.5564E-9 /
C
C     BORON
C
      DATA EB / 8.3, 25.15, 37.93, 259.37, 340.22 /
      DATA AB / 5.8365E-8, 1.0047E-8, -3.6230E-8,
     1         -7.3448E-9, 1.0220E-8, 1.6951E-9,
     2          2.0590E-8, -9.8899E-9, -6.0949E-9,
     2          2.7762E-9, 1.6499E-9, -6.7692E-10,
     3          5.7023E-9, -4.6578E-9, -9.8688E-12,
     3          1.0242E-9, -3.8855E-11, -1.7281E-10,
     4          7.3539E-10, -2.9498E-10, -3.5672E-10,
     4          1.5372E-10, 1.0036E-10, -4.8602E-11,
     5          2.5458E-10, -1.8930E-10, -2.0169E-11,
     5          4.5657E-11, 4.1076E-12, -8.7867E-12 /
      DATA AAB / 3.0952E-7, 4.7980E-8, 9.1493E-9, 1.2780E-9,
     &             4.2697E-10 /
      DATA BTB / -4.9240E-7, 1.3750E-6, -2.5387E-6,
     2           -4.1361E-8, 5.5259E-8, 9.9841E-8,
     3           -5.2910E-9, 8.6226E-9, 8.0067E-9,
     4           -1.1156E-10, -4.4920E-9, 2.3885E-8,
     5           -2.3600E-10, 1.9159E-10, 1.3125E-9 /
C
C     CARBON
C
      DATA EC / 11.26, 24.38, 47.89, 64.49, 392.08, 489.98 /
      DATA AC / 5.9848E-8, 1.1903E-8, -3.0140E-8,
     1          -1.3693E-8, 8.3748E-9, 4.0150E-9,
     2          2.8395E-8, -1.6698E-8, -2.3557E-9,
     2          3.2161E-10, 9.6016E-10, 5.2713E-10,
     3          9.0555E-9, -6.3206E-9, -1.3256E-9,
     3          1.7441E-9, 3.2680E-10, -3.8303E-10,
     4          2.7464E-9, -2.0070E-9, -2.3595E-11,
     4          4.2011E-10, -8.1600E-11, -3.9729E-11,
     5          3.9495E-10, -1.5842E-10, -1.9158E-10,
     5          8.2555E-11, 5.3899E-11, -2.6102E-11,
     6          1.4715E-10, -1.0941E-10, -1.1657E-11,
     6          2.6389E-11, 2.3742E-12, -5.0786E-12 /
      DATA AAC / 3.7442E-7, 6.0150E-8, 1.4501E-8,
     &           6.0330E-9, 6.8634E-10, 2.4679E-10 /
      DATA BTC / -6.5826E-7, 2.0521E-6, -4.4694E-6,
     2           -4.0217E-8, -2.7908E-8, 5.5499E-7,
     3           -5.3596E-9, -1.0473E-8, 1.0629E-7,
     4           -5.7710E-9, 9.9302E-9, 7.4462E-9,
     5           -5.9916E-11, -2.4124E-9, 1.2828E-8,
     6           -1.3640E-10, 1.1074E-10, 7.5862E-10 /
C
C     NITROGEN
C
      DATA EN / 14.53, 29.60, 47.45, 77.47, 97.89, 552.06, 667.03 /
      DATA AN / 4.6209E-8, 9.2264E-9, -1.2092E-8,
     1          -2.4852E-8, 5.1361E-9, 8.3068E-9,
     2          2.4369E-8, -2.2155E-9, -1.4805E-8,
     2          -4.4218E-10, 4.5211E-9, 1.7874E-10,
     3          1.2964E-8, -8.3408E-9, -2.3684E-9,
     3          2.2485E-9, 2.6234E-10, -2.6333E-10,
     4          4.6322E-9, -3.4645E-9, -3.1014E-10,
     4          8.0576E-10, 5.7791E-11, -1.4907E-10,
     5          1.5862E-9, -9.8633E-10, 7.5130E-11,
     5          3.1005E-11, -8.1970E-12, 2.6759E-11,
     6          2.3635E-10, -9.4805E-11, -1.1465E-10,
     6          4.9404E-11, 3.2255E-11, -1.5620E-11,
     7          9.2653E-11, -6.8892E-11, -7.3402E-12,
     7          1.6616E-11, 1.4949E-12, -3.1978E-12 /
      DATA AAN / 2.7367E-7, 4.4690E-8, 1.0243E-8,
     &           7.9743E-9, 5.7812E-9, 4.1073E-10,
     &           1.5539E-10 /
      DATA BTN / -4.2976E-7, 9.8352E-7, -9.5745E-7,
     2           3.0430E-8, -5.2696E-7, 2.4863E-6,
     3           3.4218E-8, -2.9155E-7, 1.2324E-6,
     4           -4.9184E-9, 6.3763E-9, 1.4917E-8,
     5           -8.5163E-9, 1.8527E-8, -7.8181E-9,
     6           -3.5856E-11, -1.4437E-9, 7.6765E-9,
     7           -8.5888E-11, 6.9728E-11, 4.7767E-10 /
C
C     OXYGEN
C
      DATA EO / 13.62, 35.12, 54.93, 77.41, 113.90, 138.12, 739.32,
     &          871.39 /
      DATA AO / 3.3559E-8, 1.3449E-8, -6.7112E-9,
     1          -1.9976E-8, 1.6214E-9, 6.5852E-9,
     2          2.4476E-8, -5.3141E-9, -7.3316E-9,
     2          -4.4515E-9, 2.4257E-9, 1.9791E-9,
     3          1.4741E-8, -8.7905E-9, -8.6099E-10,
     3          -2.4143E-10, 1.2598E-10, 6.4901E-10,
     4          6.2130E-9, -2.5047E-9, -3.0813E-9,
     4          1.3559E-9, 8.6816E-10, -4.3189E-10,
     5           2.6145E-9, -2.0276E-9, -1.6569E-10,
     5          5.0245E-10, 3.0067E-11, -9.8231E-11,
     6          1.0099E-9, -6.5165E-10, 2.8863E-12,
     6          3.0336E-11, -1.4065E-11, 4.5106E-12,
     7          1.5258E-10, -6.1203E-11, -7.4014E-11,
     7          3.1894E-11, 2.0823E-11, -1.0084E-11,
     8          6.2090E-11, -4.6167E-11, -4.9189E-12,
     8          1.1135E-11, 1.0018E-12, -2.1430E-12 /
      DATA AAO / 3.2721E-7, 4.9003E-8, 1.7523E-8,
     &           1.0270E-8, 4.0023E-9, 2.6228E-9,
     &           2.6516E-10, 1.0413E-10 /
      DATA BTO / -6.7484E-7, 2.3938E-6, -6.04388E-6,
     2           2.1747E-8, -5.4612E-7, 2.7213E-6,
     3           2.8129E-8, -3.1692E-7, 1.4381E-6,
     4           4.8134E-10, -4.3936E-8, 2.1924E-7,
     5           -1.5957E-9, -7.2485E-10, 1.9724E-8,
     6           -2.6662E-9, 2.8968E-9, 1.7246E-8,
     7           -2.3148E-11, -9.3201E-10, 4.9557E-9,
     8           -5.7557E-11, 4.6727E-11, 3.2010E-10 /
C
C----------------------------------------------------------------------
C
      INSPC = KNUCZ + 1
C
      DO 3 JSPC = 1, INSPC
      PBFI(JSPC) = 0.0
      PBFIL(JSPC) = 0.0
    3 CONTINUE
C
      KFLAG = 0
C
C     ERROR EXIT IF ELEMENT NOT AVAILABLE
C
      IF(KNUCZ .GE. 1 .AND. KNUCZ .LE. 8) GO TO 5
      KFLAG = 1
      RETURN
C
C     BRANCH ON ELEMENT TO CALL FIT EVALUATION ROUTINE WITH PROPER
C     CONSTANTS
C
    5 GO TO (10, 20, 30, 40, 50, 60, 70, 80) KNUCZ
C
   10 CALL ADBFIT(KNUCZ, PTE, PNE, EH, AH, AAH, BTH, PBFI, PBFIL)
      RETURN
   20 CALL ADBFIT(KNUCZ, PTE, PNE, EHE, AHE, AAHE, BTHE, PBFI, PBFIL)
      RETURN
   30 CALL ADBFIT(KNUCZ, PTE, PNE, ELI, ALI, AALI, BTLI, PBFI, PBFIL)
      RETURN
   40 CALL ADBFIT(KNUCZ, PTE, PNE, EBE, ABE, AABE, BTBE, PBFI, PBFIL)
      RETURN
   50 CALL ADBFIT(KNUCZ, PTE, PNE, EB, AB, AAB, BTB, PBFI, PBFIL)
      RETURN
   60 CALL ADBFIT(KNUCZ, PTE, PNE, EC, AC, AAC, BTC, PBFI, PBFIL)
      RETURN
   70 CALL ADBFIT(KNUCZ, PTE, PNE, EN, AN, AAN, BTN, PBFI, PBFIL)
      RETURN
   80 CALL ADBFIT(KNUCZ, PTE, PNE, EO, AO, AAO, BTO, PBFI, PBFIL)
      RETURN
      END
C
C#####################################################################
C
      SUBROUTINE ADBFIT(KNUCZ, PTE, PNE, PE, PA, PAA, PBT, PBFI, PBFIL)
C
C     EVALUATE THE BELFAST GROUP RATE COEFFICIENT FITTING FORMULAS.
C     NOTE THAT 0 IS RETURNED FOR TE < (IONIZATION POTENTIAL / 10).
C
      DIMENSION PE(KNUCZ), PA(6,KNUCZ), PAA(KNUCZ), PBT(3,KNUCZ)
      DIMENSION PBFI(1), PBFIL(1)
C
      ZJLEV = 1.6021E-19
C
C     CONVERT TE TO EV AND ZERO RETURN ARRAYS
C
      ZTE = 1000.0 * PTE
C
      INSPC = KNUCZ + 1
      DO 10 JSPC = 1, INSPC
      PBFI(JSPC) = 0.0
      PBFIL(JSPC) = 0.0
   10 CONTINUE
C
C     LOOP OVER ALL SPECIES EXCEPT FULLY STRIPPED
C
      DO 100 JSPC = 1, KNUCZ
C
C     RATE IS ZERO FOR KT < (IP/10)
C
      IF(ZTE .LT. PE(JSPC)/10.0) GO TO 100
C
C     BRANCH IF HIGH TEMPERATURE FIT IS REQUIRED
C
      IF(ZTE .GT. 10.0 * PE(JSPC)) GO TO 80
C
C     (IP/10) < TE < (10*IP) FITTING FORM
C
      ZX = PE(JSPC) / ZTE
      ZLOG = ALOG10(1.0 / ZX)
C
      ZSUM = PA(1,JSPC)
      DO 20 J = 2, 6
      ZSUM = ZSUM + PA(J,JSPC) * (ZLOG**(J-1))
   20 CONTINUE
C
      ZIRC = EXP(-ZX) * SQRT(1.0 / ZX) * ZSUM
      PBFI(JSPC) = PNE * ZIRC
      PBFIL(JSPC) = ZJLEV * PE(JSPC) * PBFI(JSPC)
      GO TO 100
C
C     HIGH TEMPERATURE FIT FOR TE > 10*IP
C
   80 ZX = PE(JSPC) / ZTE
C
      ZIRC = SQRT(ZX) * (PAA(JSPC) * ALOG(1.0/ZX)
     &     + PBT(1,JSPC) + PBT(2,JSPC) * ZX + PBT(3,JSPC) * ZX * ZX)
C
      PBFI(JSPC) = PNE * ZIRC
      PBFIL(JSPC) = ZJLEV * PE(JSPC) * PBFI(JSPC)
C
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADYGR(KNUCZ, PTE, PNE, PYGRI, PYGRL, KFLAG)
C
C     COMPUTE DIRECT ELECTRON IONIZATION RATES FOR HYDROGEN-LIKE
C     THROUGH NEON-LIKE CHARGE STATES OF CERTAIN HEAVY ELEMENTS
C     (PRESENTLY SCANDIUM AND IRON).  USES FORMULAE DUE TO
C     S. YOUNGER, NBS.  (JQSRT 29, 61 (1983) AND PRIVATE
C     COMMUNICATIONS).
C
C     KNUCZ= ELEMENT (Z=21 OR 26 FOR NOW)
C     PTE = ELECTRON TEMPERATURE(KEV)
C     PNE = ELECTRON DENSITY (CM-3)
C
C     PYGRI = ARRAY OF IONIZATION RATES (SEC-1)
C             (ONLY THE 10 STATES FROM H-LIKE TO NE-LIKE WILL BE
C              OVERWRITTEN BY THE SUBROUTINE)
C     PYGRL = ARRAY OF IONIZATION ENERGY LOSS RATES (WATTS/ION)
C     KFLAG = 0 FOR O.K.;
C           = 1 FOR INVALID ELEMENT (NO VALUES ARE RETURNED)
C
      COMMON / COMYGR / EION(10,3), EIONFE(10,3), EIONSC(10,3),
     &          NSS(10,3), RATE(10,3),
     &          CA(3), CB(3), CC(3), CD(3)
C
      DIMENSION PYGRI(1), PYGRL(1)
C
      DATA ((NSS(JBE,JSS),JBE=1,10),JSS=1,3)/ 1,9*2,0,0,1,7*2,
     &      4*0,1,2,3,4,5,6/
C
      DATA (CA(J),J=1,3)/7.5E-14, 5.9E-14, 8.4E-14/
      DATA (CB(J), J=1,3)/ -2.605E-14, -1.635E-14, -2.68E-14/
      DATA (CC(J), J=1,3)/ 1.195E-14, 0.82E-14, 0.672E-14/
      DATA (CD(J), J=1,3)/ -6.15E-14, -3.79E-14, -5.08E-14/
C
      DATA ((EIONSC(JBE,JSS),JBE=1,10),JSS=1,3)/ 6034., 5675.,
     &     5566., 5491., 5318., 5204., 5095., 4990., 4890., 4794.,
     &     0., 0., 1288., 1213., 1135., 1068., 993.3,
     &     914.0, 846.3, 769.3, 4*0., 1090., 1004., 922.4,
     &     838.6, 761.8, 687.4 /
C
      DATA ((EIONFE(JBE,JSS),JBE=1,10),JSS=1,3)/ 9277., 8829.,
     &     8694., 8521., 8400., 8263., 8118., 7949., 7818., 7824.,
     &     0., 0., 2046., 1950., 1852., 1755., 1679., 1566., 1482.,
     &     1394., 4*0., 1789., 1678., 1583., 1463., 1366., 1266./
C
C
C
      TEV = 1000. * PTE
C
      EION(1,1) = 0.0
      DO 20 JSS = 1, 3
      DO 20 JBE = 1, 10
      IF(KNUCZ .EQ. 21) EION(JBE,JSS) = EIONSC(JBE,JSS)
      IF(KNUCZ .EQ. 26) EION(JBE,JSS) = EIONFE(JBE,JSS)
   20 CONTINUE
C
      KFLAG = 0
      IF(EION(1,1) .NE. 0.0)GO TO 50
      KFLAG = 1
      RETURN
C
   50 DO 100 JBE = 1, 10
      ICS = KNUCZ - JBE +1
      PYGRI(ICS) = 0.0
      PYGRL(ICS) = 0.0
C
      DO 90 JSS = 1, 3
      RATE(JBE,JSS) = 0.0
      IF(NSS(JBE,JSS) .EQ. 0) GO TO 90
C
      CHI = TEV / EION(JBE,JSS)
      A = CA(JSS) * FLOAT(NSS(JBE,JSS))
      B = CB(JSS) * FLOAT(NSS(JBE,JSS))
      C = CC(JSS) * FLOAT(NSS(JBE,JSS))
      D = CD(JSS) * FLOAT(NSS(JBE,JSS))
      RATE(JBE,JSS) = FCHI(CHI,A,B,C,D) * (EION(JBE,JSS)**(-1.5)) *
     &                2.2E-06 * SQRT(CHI) * EXP(-1.0/CHI)
      PYGRI(ICS) = PYGRI(ICS) + PNE * RATE(JBE,JSS)
      PYGRL(ICS) = PYGRL(ICS) + 1.6021E-19 * PNE * EION(JBE,JSS)
     &              * RATE(JBE,JSS)
   90 CONTINUE
  100 CONTINUE
C
      RETURN
      END
C
      FUNCTION FCHI(PCHI, PA, PB, PC, PD)
C
      Z1 = PA + PB * (1.0 + 1.0 / PCHI)
      Z2 = PC - (PA + 2.0*PB + PB/PCHI)/PCHI
      Z3 = PD / PCHI
      ZCHISQ = PCHI * PCHI
      ZCHICB = ZCHISQ * PCHI
      ZALPHA = (.001193 + .9764*PCHI +.6604*ZCHISQ +
     &          .02590*ZCHICB) / (1.0 + 1.488*PCHI +
     &          .2972*ZCHISQ + .004925*ZCHICB)
      ZBETA = (-.0005725 + .01345*PCHI +.8691*ZCHISQ +
     &         .03404*ZCHICB) / (1.0 + 2.197*PCHI +
     &         .2454*ZCHISQ + .002053*ZCHICB)
      FCHI = (3.0E13 / PCHI) * (Z1 + Z2 * ZALPHA + Z3 * ZBETA)
      RETURN
      END
      SUBROUTINE ADRREC(PTE, PNE)
C
C     REV: 10/10/82 FOR AD PACKAGE
C
C     REV: 9/9/81, 6/15/78
C
C     COMPUTE RADIATIVE RECOMBINATION RATES RRAREC(JQ) FOR A GIVEN
C     ELECTRON TEMPERATURE PTE (KEV), WHERE JQ = (IONIC CHARGE + 1).
C     THESE RATES HAVE THE UNITS (SEC-1).
C     ALSO COMPUTED ARE THE RADIATION RATES ASSOCIATED WITH
C     THIS PROCESS, RADRRC(JQ).  THESE RATES ARE IN JOULES/SEC.
C
C     NOTES:
C
C     1) ELECTRONS ARE ASSUMED TO FILL IONIC SHELLS IN SIMPLE ASCENDING
C        ORDER.  THIS IS ASSUMED BOTH EXPLICITLY (THE APN ARRAY)
C        AS WELL AS IMPLICITLY THROUGH THE USE OF NVALNC VALUES
C        IN LOOP INDICES, ETC.
C     2) THROUGH USE OF THE FUNCTION EXPUND, EXP(X) IS SET EQUAL TO
C        ZERO IN THE RATE EQUATIONS FOR VALUES OF X < -45.
C     3) RATES FOR THE SPECIES JQ REFER TO TRANSFORMATIONS OF
C        THAT SPECIES INTO ADJACENT IONIZATION STATES, NOT THE
C        RATES AT WHICH THAT SPECIES IS FORMED.
C     4) UPON EVALUATION OF THE FORMULAE IT WILL BE FOUND THAT THE
C        FORMALISM USED IN XSNQ/LOKKE AND GRASBERGER (REFS 6,2)
C        IN FACT LEADS TO THE SAME EQUATIONS GIVEN IN
C        PPPL-1352/SEATON (REFS 1,3) IF THE SHIELDED CHARGE FOR
C        LEVEL N IS USED FOR 'Q'.
C
C
C     PRINCIPAL PHYSICS REFERENCES:
C
C     1) POST, JENSEN, TARTER, GRASBERGER, LOKKE, 'STEADY STATE
C        RADIATIVE COOLING RATES FOR LOW-DENSITY HIGH-TEMPERATURE
C        PLASMAS', CORRECTED PPPL-1352 (SUB. TO N.F.)
C     2) LOKKE AND GRASBERGER, 'XSNQ-U A NON-LTE EMISSION AND
C        ABSORPTION COEFFICIENT SUBROUTINE', LLL REPORT UCRL-52276
C     3) SEATON, 'RADIATIVE RECOMBINATION OF HYDROGENIC IONS',
C        MNRAS:119,81 (1959)
C     4) COX AND TUCKER, 'IONIZATION EQUILIBRIUM AND RADIATIVE
C        COOLING OF A LOW-DENSITY PLASMA', AP.J.:157,1157 (1969)
C     5) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #1, 2/11/76
C     6) XSNQ FORTRAN LISTING
C
C
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
      ZJLKEV = 1.6021E-16
C
C
C
      ZTE = PTE
C
      RRAREC(1) = 0.
      RADRRC(1) = 0.
C
C     ********** LOOP FOR ALL NON-NEUTRAL SPECIES **********
C
      DO 100 JQ = 2 , NSPC
C
      IVALNC = NVALNC(JQ)
      IVJQM1 = NVALNC(JQ - 1)
C
C
C
C     CALCULATE RADIATIVE RECOMBINATION RATES PER SHELL AND SUM THEM
C     (SUMS (FRACTIONAL VACANCY OF SHELL) * (ALPHAN) ; ALPHAN
C     FROM CORRECTED PPPL-1352).  ALSO SUM CONTINUUM AND LINE
C     RADIATION RATES DUE TO RADIATIVE RECOMBINATIONS TO EACH SHELL.
C
C
C
      ZRSUM = 0.
      ZRSUMC = 0.
      ZRSUML = 0.
C
      DO 31 JN = IVALNC,10
C
C     SKIP VALANCE SHELL IF FULL TO AVOID UNDEFINED EN, QN
C
      IF(AQN(JQ,JN) .EQ. 0.0) GO TO 31
C
      ZXN = EN(JQ,JN)/ZTE
C
      ZRTEMP = AQN(JQ,JN) * QN(JQ,JN) * (ZXN**1.5) * EXPE1(ZXN)
C
C     HIGH-LYING LEVELS CORRECTION FOR NLVLMX=10 IS 5.13
C
      IF(JN .EQ. 10) ZRTEMP = ZRTEMP * 5.13
      ZRSUM = ZRSUM + ZRTEMP
C
      ZRTMPL = ENM(JQ-1,IVJQM1,JN) * ZRTEMP
      ZRSUML = ZRSUML + ZRTMPL
C
      ZRTMPC = AQN(JQ,JN) * QN(JQ,JN) * SQRT(ZXN) * EN(JQ,JN)
      ZRSUMC = ZRSUMC + ZRTMPC
   31 CONTINUE
C
      RRAREC(JQ) = PNE * 5.20E-14 * ZRSUM
      IF(ZRSUMC .LT. 1.E-07) ZRSUMC = 0.
      ZRSUMC = 8.32E-30 * ZRSUMC
C
      IF(ZRSUML .LT. 1.E-7) ZRSUML = 0.
      ZRSUML = ZJLKEV * 5.20E-14 * ZRSUML
      RADRRC(JQ) = PNE * (ZRSUMC + ZRSUML)
C
C
  100 CONTINUE
C
C
      RETURN
      END
      SUBROUTINE ADDREC(PTE, PNE)
C7600      LCM DIEREC
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 11/11/81 ADD CDNN AND CDNM MULTIPLIERS
C
C     REV: 10/23/78 (REMOVE 'C' FROM ZDSUM0 DEFINING STATEMENT TO FIX
C                    N-N DIELECTRONIC RATE BUG)
C     REV: 6/15/78
C
C     COMPUTES DIELECTRONIC RECOMBINATION RATE RDIREC(JQ)
C     FOR A GIVEN ELECTRON TEMPERATURE PTE (KEV) AND ELECTRON
C     DENSITY PNE (CM-3), WHERE JQ = (IONIC CHARGE +1).
C     THESE RATES HAVE THE UNITS (SEC-1).  SINCE
C     THE RATE COEFF DEPEND ON NE, THE TOTAL RECOMBINATION RATE
C     PER ION IS NOT STRICTLY A LINEAR FUNCTION OF NE.  HOWEVER,
C     THE DIELECTRONIC COLLISIONAL-INTERRUPTION CORRECTION IS
C     ONLY WEAKLY DEPENDENT ON THE DENSITY, AND THUS THE NON-LINEARITY
C     IS GENERALLY SMALL.
C
C     ALSO COMPUTED ARE THE RADIATION RATES ASSOCIATED WITH THIS
C     PROCESS, RADDRC(JQ), WHICH INCLUDE CONTRIBUTIONS FROM
C     THE RADIATIVE DECAYS OF BOTH THE EXCITED AND CAPTURED ELECTRONS.
C     THESE RATES ARE IN JOULES/SEC.
C
C     NOTES:
C
C     1) ELECTRONS ARE ASSUMED TO FILL IONIC SHELLS IN SIMPLE ASCENDING
C        ORDER.  THIS IS ASSUMED BOTH EXPLICITLY (THE APN ARRAY)
C        AS WELL AS IMPLICITLY THROUGH THE USE OF NVALNC VALUES
C        IN LOOP INDICES, ETC.
C     2) THROUGH USE OF THE FUNCTION EXPUND, EXP(X) IS SET EQUAL TO
C        ZERO IN THE RATE EQUATIONS FOR VALUES OF X < -45.
C     3) RATES FOR THE SPECIES JQ REFER TO TRANSFORMATIONS OF
C        THAT SPECIES INTO ADJACENT IONIZATION STATES, NOT THE
C        RATES AT WHICH THAT SPECIES IS FORMED.
C
C
C     PRINCIPAL PHYSICS REFERENCES:
C
C     1) POST, JENSEN, TARTER, GRASBERGER, LOKKE, 'STEADY STATE
C        RADIATIVE COOLING RATES FOR LOW-DENSITY HIGH-TEMPERATURE
C        PLASMAS', CORRECTED PPPL-1352 (SUB. TO N.F.)
C     2) BURGESS, 'A GENERAL FORMULA FOR THE ESTIMATION OF
C        DIELECTRONIC RECOMBINATION COEFFICIENTS IN LOW-DENSITY
C        PLASMAS', AP.J. 141, 1588 (1965).
C     3) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #1, 2/11/76
C     4) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #4, 7/29/76
C     5) GRASBERGER, XSNQ - TOKAMAK MEMORANDUM #10, 1/5/77
C     6) MERTS, COWAN, MAGEE, 'THE CALCULATED POWER OUTPUT
C        FROM A THIN IRON-SEEDED PLASMA', LOS ALAMOS REPORT LA-6220-MS
C     7) XSNQ FORTRAN LISTING
C     8) ANSARI, ELWERT AND MUCKLICH, 'ON DIELECTRONIC RECOMBINATION',
C        ZEITSCHRIFT NATURFORSCHUNG 25A, 1781 (1970).
C
C
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
      COMMON /ADDIEL/ CDNN(100), CDNM(100), LDRMLT, YHNNA, YHNNB,
     &                YHNNC, YHNMA, YHNMB, YHNMC

C
      ZJLKEV = 1.6021E-16
C
C     CALL ADDRCC TO FILL CDNN AND CDNM MULTIPICATIVE CONSTANT ARRAYS
C
      CALL ADDRCC
C
      ZTE = PTE
      ZNE = PNE
C
      ZNUC = FLOAT(NUCZ)
C
      RDIREC(1) = 0.
      RADDRC(1) = 0.
      RDIREC(NSPC) = 0.
      RADDRC(NSPC) = 0.
C
      IF(NUCZ .EQ. 1) RETURN
C
C
C     ********** LOOP FOR ALL SPECIES EXCEPT NEUTRAL **********
C              ********** AND FULLY IONIZED **********
C
      DO 100 JQ = 2 , NUCZ
C
C
      ZCHG = FLOAT(JQ-1)
      IVALNC = NVALNC(JQ)
C
C
C
C     ************* DIELECTRONIC RECOMBINATION RATE *************
C
C     (EQUATIONS AND VARIABLES FOLLOW PPPL-1352/BURGESS FORMALISM)
C
C
C
      ZQ = ZCHG
      ZQP1 = ZCHG + 1
C
      ZBQ = SQRT(  (ZQ * (ZQ+1.)**5) / (ZQ*ZQ + 13.4)  )
      ZNT = (4.77E+18 * (ZQ**6) * SQRT(ZTE) / ZNE)**0.142857143
      ZDNN = .0068 * ZNT * (1.+.001*ZNT) / (1.+.00002*ZNT*ZNT)
      ZZZ = .0015 * ( (ZQ+1.) * ZNT )**2
      ZDNM = ZZZ / (1. + ZZZ)
      ZZZNN = 7.59E-14 * ZBQ * ZDNN / (ZTE**1.5)
      ZZZNM = 7.59E-14 * ZBQ * ZDNM / (ZTE**1.5)
      ZZZA = 1. + .015 * (ZQ**3) / ( (ZQ+1.)**2 )
      ZZZET = ( (ZQ + 1.)**2 ) / (73.5 * ZZZA * ZTE)
C
C
C     COMPUTE N - N CONTRIBUTION TO SUM FIRST
C
C
      ZEIJNN = ENN(JQ) / (.0136 * ZQP1 * ZQP1)
      ZYNN = (ZQ + 1.) * ZEIJNN
      ZAYNN = SQRT(ZYNN) / (1. + .105 * ZYNN + .015 * ZYNN**2)
C
      ZDSUM0 = FNN(JQ) * ZAYNN * EXPUND(-ZZZET*ZEIJNN , -45.)
C
C
C     SUM OVER N - M TRANSITIONS (NOTE THAT TRANSITIONS FROM
C     INNER SHELLS TO THE VALENCE SHELL AND ABOVE ARE INCLUDED)
C
C
      ZDSUM = 0.
      ZDSUMR = 0.
C
      DO 40 JN = 1 , IVALNC
      IIJM = IVALNC
      IF(JN .EQ. IVALNC .OR.
     &          AQN(JQ,IVALNC) .EQ. 0.) IIJM = IVALNC + 1
C
      DO 40 JM = IIJM , 10
C
      ZEIJ = ENM(JQ,JN,JM) / (.0136 * ZQP1**2)
      ZY = (ZQ + 1.) * ZEIJ
      ZAY = SQRT(ZY) / (2. + .420 * ZY + .06 * ZY**2)
C
      ZDTEMP = FNM(JQ,JN,JM) * ZAY * EXPUND(-ZZZET*ZEIJ , -45.)
C
      ZDSUM = ZDSUM + ZDTEMP
      ZDSUMR = ZDSUMR + ENM(JQ,JN,JM) * ZDTEMP
   40 CONTINUE
C
      ZDNNRT = ZZZNN * ZDSUM0
      ZDNMRT = ZZZNM * ZDSUM
CW    WRITE(59,909)JQ,ZDNNRT,ZDNMRT
CW  909 FORMAT(1X,'JQ =',I3,3X,'ZDNNRT =',1PE12.2,'  ZDNMRT =',E12.2)
      RDIREC(JQ) = PNE * (CDNN(JQ) * ZDNNRT + CDNM(JQ) * ZDNMRT)
C
      ITEMP = NVALNC(JQ-1)
      ZIQ = EN(JQ,ITEMP)
      ZTEMP = CDNN(JQ) * (ZIQ + ENN(JQ)) * ZDNNRT +
     &              CDNM(JQ) * (ZIQ * ZDNMRT + ZZZNM * ZDSUMR)
      IF(ZTEMP .LT. 1.E-21) ZTEMP = 0.
      RADDRC(JQ) = PNE * ZJLKEV * ZTEMP
C
C
C
  100 CONTINUE
C
C
      RETURN
      END
      SUBROUTINE ADDRCC
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     SETUP THE DIELECTRONIC RECOMBINATION RATE COEFFICIENT
C     MULTIPLICATIVE CONSTANT ARRAYS CDNN AND CDNM ACCORDING
C     TO THE LDRMLT SWITCH.
C
C     LDRMLT = 0     USE CDNN, CDNM ARRAYS AS GIVEN
C     LDRMLT = 1     CDNN, CDNM ARRAYS SET TO 1.0
C     LDRMLT = 2     USE Y. HAHN FACTORS TO SETUP CDNN, CDNM
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADDIEL / CDNN(100), CDNM(100), LDRMLT,
     &                  YHNNA, YHNNB, YHNNC, YHNMA, YHNMB, YHNMC
C
C
      IF(LDRMLT .EQ. 0) RETURN
C
      DO 20 JSPC = 1, 100
      CDNN(JSPC) = 1.0
      CDNM(JSPC) = 1.0
   20 CONTINUE
C
      IF(LDRMLT .EQ. 1) RETURN
C
C
      DO 100 JSPC = 2, NUCZ
C
      ZNELEC = FLOAT(NUCZ - JSPC + 1)
C
      CDNN(JSPC) = YHNNA * EXP(-YHNNB*SQRT(ZNELEC) ) *
     &                  (FLOAT(NUCZ)**YHNNC)
C
      CDNM(JSPC) = YHNMA * EXP(-YHNMB*SQRT(ZNELEC)) *
     &                  (FLOAT(NUCZ)**YHNMC)
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADECEX(PTE, PNE)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     COMPUTE RADIATION RATES RADCLX(JQ) (JOULE/SEC) DUE TO
C     COLLISIONAL EXCITATION FOR EACH IONIC SPECIES JQ = (IONIC
C     CHARGE + 1).  ALL EXCITATIONS ARE ASSUMED TO RESULT IN
C     RADIATIVE DECAYS.
C
C     REFERENCES:
C     1) POST ET AL. PPPL-1352 (CORRECTED)
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
      COMMON / CEXDAT / ZDNN(28), ZGAMNN(28)
C
C
      DATA (ZDNN(J),J=1,28) / 0., 0.,
     & 1.2, .91, .89, .92, .86, .87, .85, 0.,
     & .72, .78, .76, .78, .75, .70, .68, .65, .65,
     & 8 * .70, 0. /
C
      DATA (ZGAMNN(J),J=1,28) / 0., 0.,
     & .54, .77, .58, .57, .41, .63, .66, 0.,
     & .97, .96, .88, .86, .87, .85, .83, .81,
     & 9 * .80, 0. /
C
      ZJLKEV = 1.6021E-16
C
C
C
      ZTE = PTE
      ZNUC = FLOAT(NUCZ)
C
C     ******** LOOP FOR ALL SPECIES EXCEPT FULLY IONIZED ********     
C
      DO 100 JQ = 1 , NUCZ
C
      ZCHG = FLOAT(JQ - 1)
      IVALNC = NVALNC(JQ)
      IBOUND = NUCZ - JQ + 1
C
C     ***** N - N RATE *****
C
      ZZNN = 0.
      IF(FNN(JQ) .EQ. 0.) GO TO 30
C
      ZXNN = ENN(JQ) / ZTE
C
      IF(JQ .NE. 1) GO TO 20
      ZTEMP = .06 * (SQRT(ZXNN) - 2.) / (1. + ZXNN)
      GO TO 22
C
   20 IF(IBOUND .LE. 28) GO TO 21
      ZTEMP = .7 * (1. - .8 / ZCHG)
      GO TO 22
C
   21 ZTEMP = ZDNN(IBOUND) * (1. - ZGAMNN(IBOUND) / ZCHG)
C
   22 ZGNN = ZTEMP + .276 * EXPE1(ZXNN)
C
      ZZNN = FNN(JQ) * EXPUND(-ZXNN , -45.) * ZGNN
C
C     ***** SUM N - M RATES FROM ALL OCCUPIED SHELLS *****
C
   30 ZNMSUM = 0.
C
C
      DO 31 JN = 1 , IVALNC
C
      ZJN = FLOAT (JN)
      IIJM = IVALNC
      IF(JN .EQ. IVALNC .OR.
     &        AQN(JQ,IVALNC) .EQ. 0.) IIJM = IVALNC + 1
C
      DO 31 JM = IIJM , 10
      ZJM = FLOAT(JM)
C
      ZXNM = ENM(JQ,JN,JM) / ZTE
      ZTEMP = ZJM * (ZJM - ZJN) * (1. + (1. - 2./ZNUC) * ZXNM) / 20.
      ZGNM = .19 * (1. + .9 * (1. + ZTEMP) * EXPE1(ZXNM))
C
C     NOTE THAT FNM HAS APN AND AQN FACTORS ALREADY INCLUDED
C
      ZZNM = FNM(JQ,JN,JM) * EXPUND(-ZXNM , -45.) * ZGNM
      ZNMSUM = ZNMSUM + ZZNM
   31 CONTINUE
C
C
C
      ZTEMP = (ZZNN + ZNMSUM) / SQRT(ZTE)
      IF(ZTEMP .LT. 1.E-10) ZTEMP = 0.
      RADCLX(JQ) = PNE * ZJLKEV * 4.99E-10 * ZTEMP
C
  100 CONTINUE
C
C     ***************END OF JQ SPECIES LOOP ***************
C
      RADCLX(NSPC) = 0.
C
      RETURN
      END
      SUBROUTINE ADBREM(PTE, PNE)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 6/15/78
C
C     INPUT : PTE, NUCZ, NSPC
C     OUTPUT : RADBRM
C
C     COMPUTE BREMSSTRAHLUNG RADIATION RADBRM(JQ) (JOULE/SEC)
C     FOR EACH IONIC SPECIES JQ = (IONIC CHARGE + 1).
C
C     REFERENCES:
C     1) POST ET AL. PPPL-1352 (CORRECTED)
C     2) KARZAS AND LATTER
C     3) TUCKER
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10), FNN(100), ENN(100)
C
      COMMON / ADELEC / RCLION(100), RRAREC(100), RDIREC(100),
     &                  CIZLOS(100), RADRRC(100), RADDRC(100),
     &                  RADCLX(100), RADBRM(100)
C
C
      ZTE = PTE
C
C     USE FIXED FREE-FREE GAUNT FACTOR ZGFF = 1.2
C
      ZGFF = 1.2
C
      DO 100 JQ = 1 , NSPC
      ZCHG = FLOAT(JQ - 1)
      RADBRM(JQ) = PNE * 4.85E-31 * ZCHG * ZCHG * SQRT(ZTE) * ZGFF
  100 CONTINUE
C
      RETURN
      END
      FUNCTION EXPE1(PX)
C
C     REV: 6/15/78
C
C     CALCULATE EXP(PX) * E1(PX), WHERE THE EXPONENTIAL INTEGRAL
C     E1(X) = INTEGRAL FROM X TO INFINITY OF (EXP(-T)/T)DT.
C
C     REFERENCE: 'HANDBOOK OF MATHEMATICAL FUNCTIONS', ABRAMOWITZ
C     AND STEGUN, NATIONAL BUREAU OF STANDARDS
C
      IF(PX.GE.1.)GO TO 2
      IF(PX.GT.1.E-4)GO TO 1
C
      EXPE1 = EXP(PX) * ( - ALOG(PX) - .57721566 + PX )
      RETURN
C
C
    1 EXPE1 = EXP(PX) * ( - ALOG(PX) - .57721566 + .99999193*PX
     &         - .24991055*PX**2 + .05519968*PX**3 - .00976004*PX**4
     &         + .00107857*PX**5 )
C
      RETURN
C
C
    2 EXPE1 = (PX + 2.334733 + .250621/PX) /
     &         (PX**2 + 3.330657*PX + 1.681534)
C
      RETURN
      END
C
C
C********************************************************************
C
C
      FUNCTION EXPUND(PX,PCHK)
C
C     REV: 6/15/78
C
C     PREVENT UNDERFLOWS IN RATE CALCULATIONS BY RETURNING
C     0. FOR EXP(PX) IF PX < PCHK.
C
      IF( PX - PCHK ) 1 , 1 , 2
    1 EXPUND = 0.
      RETURN
C
    2 EXPUND = EXP(PX)
      RETURN
C
      END
      SUBROUTINE ADBCXR(PNNEUT, PVNEUT, KNCXB, KCXOPT, KERR, KVUNIT)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     REV: 7/24/79; ZERO RADBCX BEFORE KNCXB=0 TEST & RETURN
C     REV: 5/10/79; ADD RADIATION DUE TO BEAM CHARGE EXCHANGE VIA
C                   RADBCX AND BEQRAD
C
C     COMPUTE RECOMBINATION RATE DUE TO CHARGE EXCHANGE WITH NEUTRAL
C     BEAM IONS.
C
C          NPSC = NUMBER OF IMPURITY SPECIES (NUCLEAR Z + 1)
C
C          PNNEUT = ARRAY OF ANGULAR AVERAGED NEUTRAL DENSITIES FOR
C                   EACH BEAM ENERGY COMPONENT
C          PVNEUT = ARRAY OF VELOCITIES FOR EACH BEAM COMPONENT (CM/S)
C          KNCXB = NUMBER OF ENERGY COMPONENTS IN THE BEAM
C                 (KNCXB = 0 TURNS OFF RCXREC)
C          KCXOPT = 1, 2, OR 3 TO SELECT CROSS-SECTIONS TO BE USED
C          RCXREC = ARRAY OF CHARGE EXCHANGE RECOMBINATION RATES FOR
C                   EACH SPECIES, SEC-1 UNITS (E.G. RATE PER ION)
C          KERR = ERROR FLAG FROM ADCXXS ROUTINE (0 = O.K.)
C          RADBCX = RADIATION RATE DUE TO CASCADE AFTER
C                   CHARGE EXCHANGE TO EXCITED LEVEL (WATTS / ION)
C
C
      COMMON / ADSDAT / NUCZ, NSPC, APN(100,10), AQN(100,10),
     &                  NVALNC(100), SIGMA(10,10), QN(100,10),
     &                  EN(100,10), EIN(100,5), ENM(100,5,10),
     &                  FNM(100,5,10),FNN(100),ENN(100)
C
      COMMON / ADNEUT / RCXREC(100), RADBCX(100)
C
      DIMENSION PNNEUT(1), PVNEUT(1)
      DIMENSION ZSIGMA(100)
C
      ZJLKEV = 1.602E-16
C
      INSPC = NSPC
C
C     ZERO RATE ARRAYS FOR ACCUMULATION OF CONTRIBUTIONS FROM EACH BEAM
C     COMPONENT
C
      DO 10 JSPC = 1 , 100
      RCXREC(JSPC) = 0.
      RADBCX(JSPC) = 0.
   10 CONTINUE
C
C     KNCXB = 0 OR KCXOPT = 0 TO TURN OFF C-X RECOMBINATION
C
      IF(KNCXB .EQ. 0 .OR. KCXOPT .EQ. 0) RETURN
C
C     OBTAIN CROSS-SECTIONS FOR EACH VELOCITY (ENERGY) COMPONENT
C     IN THE BEAM AND COMPUTE TOTAL RECOMBINATION RATE.
C
      DO 60 JB = 1 , KNCXB
C
      CALL ADCXXS(PVNEUT(JB), KVUNIT, ZVCMS, ZSIGMA, INSPC,
     &                KCXOPT, KERR)
      ZNFLUX = ZVCMS * PNNEUT(JB)
C
      DO 50 JSPC = 2 , INSPC
      RCXREC(JSPC) = RCXREC(JSPC) + ZNFLUX * ZSIGMA(JSPC)
   50 CONTINUE
C
   60 CONTINUE
C
      DO 70 JSPC = 2 , INSPC
      ITEMP = NVALNC(JSPC-1)
      RADBCX(JSPC) = ZJLKEV * EN(JSPC,ITEMP) * RCXREC(JSPC)
   70 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADCXXS(PVEL, KVUNIT, PVCMS, PSIGMA, KNSPC,
     &                     KOPT, KERR)
C
C     REV: 10/10/82 FOR NEW AD PACKAGE
C
C     SUPPLIES CROSS-SECTIONS FOR CHARGE EXCHANGE BETWEEN NEUTRAL
C     HYDROGEN AND MULTI-CHARGED IONS.  THE VARIABLE KOPT ALLOWS
C     THE SELECTION OF RESULTS FROM ONE OF THE VARIOUS PUBLISHED
C     CALCULATIONS.
C
C     PVEL =  RELATIVE VELOCITY, UNITS AS SELECTED BY KVUNIT
C     KVUNIT = UNITS FLAG FOR INPUT VELOCITY
C                 1 = CM/SEC
C                 2 = KEV/AMU
C                 3 = ATOMIC UNITS
C     PSIGMA = ARRAY OF CROSS-SECTION VALUES RETURNED, CM**2.
C              PSIGMA(I) = CROSS-SECTION FOR ION WITH CHARGE = I-1,
C              FOR 2 <= I <= KNSPC
C     KNSPC = HIGHEST SPECIES NUMBER, = MAX Z + 1
C     KOPT = SWITCH FOR DIFFERENT RESULTS, AS FOLLOWS:
C            1 = OLSON & SALOP, PHYS. REV. A 14, 579 (1976)
C                (ABSORBING SPHERE, V <= 1.E08 CM/SEC)
C            2 = GROZDANOV & JANEV, PHYS. REV. A 17, 880 (1978)
C                (TUNNELING, 2.E06 <= V <= 8.E08 CM/SEC)
C            3 = OLSON & SALOP, PHYS. REV. A 16, 531 (1977)
C                (CLASSICAL, (2 - 7)E08 CM/SEC)
C     KERR = ERROR FLAG. 0 = O.K., <0 = INVALID INPUT, OTHERS IN TEXT
C
C*********************************************************************
C
      DIMENSION PSIGMA(1)
      COMMON / CXXDAT / ZFITV(4), ZFITLT(4), ZFITGT(4)
C
      DATA ZFITV/ 2.7E08, 3.0E08, 4.0E08, 5.0E08/
      DATA ZFITLT/1.220988, 1.277249, 1.868400, 2.575095/
      DATA ZFITGT/1.000000, 1.050777, 1.221593, 1.668531/
C
C---------------------------------------------------------------------
C
      KERR = 0
C
C     CONVERT INPUT VELOCITY PARAMETER
C
      IF(KVUNIT .NE. 1) GO TO 2
      ZVCMS = PVEL
      ZVKEVM = (ZVCMS / 4.39E07)**2
      ZVAU = ZVCMS / 2.19E08
      GO TO 5
C
    2 IF(KVUNIT .NE. 2) GO TO 3
      ZVKEVM = PVEL
      ZVCMS = 4.39E07 * SQRT(ZVKEVM)
      ZVAU = ZVCMS / 2.19E08
      GO TO 5
C
    3 IF(KVUNIT .NE. 3) GO TO 4
      ZVAU = PVEL
      ZVCMS = 2.19E08 * ZVAU
      ZVKEVM = (ZVCMS / 4.39E07)**2
      GO TO 5
C
    4 KERR = - 1
      RETURN
C
C     RETURN VELOCITY IN CM/S UNITS
C
    5 PVCMS = ZVCMS
C
C     NEUTRAL SPECIES CROSS-SECTION = 0
C
      PSIGMA(1) = 0.
C
      GO TO (100, 200, 300), KOPT
      KERR = -2
      RETURN
C
C--------------------------------------------------------------------
C
C     KOPT = 1 ; OLSON & SALOP LOW - V ABSORBING SPHERE
C     NOMINALLY VALID ONLY FOR V < 1.0E08 CM/SEC
C
  100 DO 110 JSPC = 3 , KNSPC
C
CW    WRITE(6,198)JSPC
CW  198 FORMAT(1X,/,'JSPC =',I3)
C
      ZCHG = FLOAT(JSPC - 1)
      ZL = 2.648 / SQRT(ZCHG)
      ZC = 2.864E-04 * ZCHG * (ZCHG - 1.) * ZVAU
      IF(ZC .LT. .541341133/(ZL**2) ) GO TO 102
      KERR = 1
      RETURN
C
C     ITERATIVE SOLUTION TO EQUATION (10)
C
C     (FIRST GUESS, FROM LINEAR FIT TO FIG. 3)
C
  102 ZR = 9. + 11. * ZCHG / 45.
      ZR = - ALOG( ZC / (ZR**2) ) / ZL
      ZFR = ZR * ZR * EXP(- ZL * ZR)
      ICOUNT = 0
C
  105 IF(ZR .LT. 2./ZL) GO TO 180
      IF(ICOUNT .GT. 100) GO TO 181
      ICOUNT = ICOUNT + 1
C
CW    WRITE(6,199)ZR
  199 FORMAT(1X,'ZR=',F10.2)
C
      ZR = ZR + (ZC/ZFR - 1.) / (2./ZR - ZL)
      ZFR = ZR * ZR * EXP (- ZL * ZR)
C
      IF( ABS(ZFR - ZC)/ZC .GE. .001) GO TO 105
C
C     ITERATION COMPLETE TO 0.1% IN VELOCITY; COMPUTE SIGMA
C
      PSIGMA(JSPC) = 3.14159 * (0.529E-08 * ZR)**2
  110 CONTINUE
C
C     FILL IN Z = 1 CROSS SECTION WITH Z = 2 VALUE
C
      PSIGMA(2) = PSIGMA(3)
      RETURN
C
  180 KERR = 2
      RETURN
  181 KERR = 3
      RETURN
C
C----------------------------------------------------------------
C
C     KOPT = 2; GROZDANOV & JANEV
C     RESULTS SHOWN ONLY FOR Z = 10, 20, AND 30, 2.E06 < V < 8.E08
C     CM/SEC; ALL OTHER POINTS BY INTERPOLATION OR EXTRAPOLATION
C
  200 DO 210 JSPC = 2 , KNSPC
      ZCHG = FLOAT(JSPC - 1)
C
C     EMPIRICAL FIT TO FIGURE 3, CENTERED AT V = 2.0E08 CM/SEC
C
      ZSLOPE = - (.25 + .04 * ZCHG) * 1.0E-14
      Z0 = .11 * ZCHG * 1.0E-14
      PSIGMA(JSPC) = AMAX1(0., Z0 + ZSLOPE * ALOG10(ZVCMS/2.E08))
C
  210 CONTINUE
      RETURN
C
C----------------------------------------------------------------------
C
C
C     KOPT = 3 ; OLSON & SALOP CLASSICAL TRAJECTORY.  USES A
C     POWER LAW FIT TO FIG. 1 AT 4 VELOCITIES IN THE RANGE
C     2.7E08 < V < 5.E08 CM/SEC.  INPUT VELOCITIES OUTSIDE THIS
C     RANGE YIELD DATA FOR THE APPROPRIATE LIMITING VELOCITY.  THE
C     LAST GIVEN DATA CURVE IS FOR Z=36, HENCE DATA FOR Z>36
C     REPRESENTS AN EXTRAPOLATION.
C
  300 CONTINUE
      Z3V = AMAX1(2.7E08, ZVCMS)
      Z3V = AMIN1(5.0E08, Z3V)
C
      IX = IFIX(Z3V/1.E08) - 1
      IF(IX .GT. 3) IX = 3
      IXP1 = IX + 1
C
C     USE LINEAR FIT TO Z=8 CURVE TO FIND X-SECTION FOR Z=8 AT GIVEN
C     VELOCITY
C
      ZSIG8 = 7.43E-15 - 1.26E-23 * Z3V
C
C     INTERPOLATE TO FIND POWER LAW EXPONENTS FOR THIS VELOCITY; ONE
C     FOR Z<8, THE OTHER FOR Z>8.
C
      ZALT = ZFITLT(IX) + (ZFITLT(IXP1) - ZFITLT(IX)) *
     &                 ( (Z3V-ZFITV(IX)) / (ZFITV(IXP1) - ZFITV(IX)) )
      ZAGT = ZFITGT(IX) + (ZFITGT(IXP1) - ZFITGT(IX)) *
     &                 ( (Z3V-ZFITV(IX)) / (ZFITV(IXP1)-ZFITV(IX)) )
C
C     FIT PARAMETERS FOUND FOR THIS VELOCITY; COMPUTE X-SECTIONS
C
      DO 315 JSPC = 2 , KNSPC
      ZDIV8 = FLOAT(JSPC - 1) / 8.
      IF(ZDIV8 .GE. 1.) GO TO 310
      PSIGMA(JSPC) = ZSIG8 * ZDIV8**ZALT
      GO TO 315
  310 PSIGMA(JSPC) = ZSIG8 * ZDIV8**ZAGT
  315 CONTINUE
C
      RETURN
C
C---------------------------------------------------------------------
C
      END