C**********************************************************************
C
C  TDFLIB - LIBRARY OF SUBPROGRAMS TO READ THERMONUCLEAR DATA FILE TDF
C
C  DEVELOPED AND CODED BY: 
C     S. I. WARSHAW
C     NUCLEAR DATA GROUP, L-298 
C     CP DIVISION, PHYSICS DEPARTMENT
C     LAWRENCE LIVERMORE NATIONAL LABORATORY
C     LIVERMORE, CA 94550 USA
C
C  VERSION DP3E,  12 NOVEMBER 1991
C
C**********************************************************************
C
C  QUESTIONS OR COMMENTS ABOUT THE DATA IN TDF OR THE
C  ASSOCIATED ROUTINES SHOULD BE ADDRESSED TO THE LLNL
C  NUCLEAR DATA GROUP AS FOLLOWS.
C
C  BY SURFACE MAIL:
C     GROUP LEADER FOR NUCLEAR DATA, L-298
C     LAWRENCE LIVERMORE NATIONAL LABORATORY
C     P.O. BOX 808
C     LIVERMORE, CA  94550 USA
C
C  BY ELECTRONIC MAIL:  NDG@LLNL.GOV
C
C  BY TELEPHONE: (510) 422-9668
C
C  BY TELEFAX: (510) 422-9523
C
C**********************************************************************
C
C  COPYRIGHT NOTICE:
C
C  COPYRIGHT 1991, THE REGENTS OF THE UNIVERSITY OF
C  CALIFORNIA.  ALL RIGHTS RESERVED.  THIS WORK WAS PRODUCED
C  AT THE UNIVERSITY OF CALIFORNIA, LAWRENCE LIVERMORE
C  NATIONAL LABORATORY (UCLLNL) UNDER CONTRACT NO.
C  W-7405-ENG-48 (CONTRACT 48) BETWEEN THE U. S. DEPARTMENT
C  OF ENERGY (DOE) AND THE REGENTS OF THE UNIVERSITY OF
C  CALIFORNIA (UNIVERSITY) FOR THE OPERATION OF UCLLNL. 
C  COPYRIGHT IS RESERVED TO THE UNIVERSITY FOR PURPOSES OF
C  CONTROLLED DISSEMINATION, COMMERCIALIZATION THROUGH
C  FORMAL LICENSING, OR OTHER DISPOSITION UNDER THE TERMS OF
C  CONTRACT 48;  DOE POLICIES, REGULATIONS AND ORDERS; AND
C  U. S. STATUTES.  THE RIGHTS OF THE FEDERAL GOVERNMENT ARE
C  RESERVED UNDER CONTRACT 48 SUBJECT TO THE RESTRICTIONS
C  AGREED UPON BY DOE AND UNIVERSITY AS ALLOWED UNDER DOE
C  ACQUISITION LETTER 88-1.
C
C
C  DISCLAIMER:
C
C  THIS COMPUTER SOFTWARE AND ANCILLARY DOCUMENTATION WAS
C  PREPARED AS AN ACCOUNT OF WORK SPONSORED BY AN AGENCY
C  OF THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED
C  STATES GOVERNMENT NOR THE UNIVERSITY OF CALIFORNIA NOR
C  ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR
C  IMPLIED, OR ASSUMES ANY LIABILITY OR RESPONSIBILITY FOR
C  THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY
C  INFORMATION, APPARATUS, PRODUCT, OR PROCESS DISCLOSED, OR
C  REPRESENTS THAT ITS SPECIFIC COMMERCIAL PRODUCTS,
C  PROCESS, OR SERVICE BY TRADE NAME, TRADEMARK,
C  MANUFACTURER, OR OTHERWISE,  DOES NOT NECESSARILY
C  CONSTITUTE OR IMPLY ITS ENDORSEMENT, RECOMMENDATION, OR
C  FAVORING BY THE UNITED STATES GOVERNMENT OR THE
C  UNIVERSITY OF CALIFORNIA. THE VIEWS AND OPINIONS OF THE
C  AUTHORS EXPRESSED HEREIN DO NOT NECESSARILY STATE OR
C  REFLECT THOSE OF THE UNITED STATES GOVERNMENT OR THE
C  UNIVERSITY OF CALIFORNIA, AND SHALL NOT BE USED FOR
C  ADVERTISING OR PRODUCT ENDORSEMENT PURPOSES.
C
C**********************************************************************
C
C  REFERENCES:
C
C  1. "A GENERALIZED ENERGY EXCHANGE KERNEL FOR INELASTIC
C     NEUTRON SCATTERING AND THERMONUCLEAR REACTIONS",
C     M. M. R. WILLIAMS, JOURNAL OF NUCLEAR ENERGY VOLUME 25,
C     PP. 489-501 (1971)
C
C  2. "FUSION NEUTRON ENERGIES AND SPECTRA", H. BRYSK, PLASMA
C     PHYSICS VOLUME 15, PP.611-617 (1973)
C
C  3. "EVALUATION OF CHARGED PARTICLE REACTIONS FOR FUSION
C     APPLICATIONS", R. M. WHITE, D. A. RESLER AND
C     S. I. WARSHAW, LAWRENCE LIVERMORE NATIONAL LABORATORY
C     REPORT UCRL-JC-107158, 7 PP. (MAY, 1991);  ALSO IN THE
C     PROCEEDINGS OF THE INTERNATIONAL CONFERENCE ON NUCLEAR
C     DATA FOR SCIENCE AND TECHNOLOGY, JULICH, GERMANY,
C     13-17 MAY 1991, TO BE PUBLISHED.
C    
C**********************************************************************
C
C  NOTES:
C
C  THE FOLLOWING FIVE ROUTINES ARE INTENDED FOR DIRECT USE BY THE
C  USER'S APPLICATIONS CODE:
C                       TDFRDR
C                       SVBLU 
C                       MXAVLU
C                       SPECLU
C                       TDMASS
C
C  THE FOLLOWING FOUR SUBORDINATE ROUTINES ARE CALLED BY THE
C  ROUTINES LISTED ABOVE:
C                       TDCLAG
C                       TDCHLU
C                       TDMCHL
C                       TDSELU
C
C  THE FOLLOWING ARE THE NAMES OF THE THREE LABELED COMMON BLOCKS
C  THAT ARE SET UP IN MEMORY:
C                       TDFCON
C                       TDFSVB
C                       TDFSPE
C
C  THESE ROUTINES ARE CODED IN ANSI STANDARD FORTRAN 77.
C  THEY ARE INTENDED TO BE RUN ON CPUS THAT USE 32-BIT WORD
C  ARCHITECTURE, HENCE ALL FLOATING POINT VARIABLES AND
C  NUMERICAL CONSTANTS HEREIN ARE TYPED DOUBLE PRECISION.
C  IF IT IS DESIRED TO USE THESE ROUTINES ON CPUS EMPLOYING
C  64-BIT WORD LENGTHS, PLEASE CONTACT THE LLNL NUCLEAR
C  DATA GROUP, AS INDICATED ABOVE.
C
C  IN ARGUMENT DESCRIPTIONS WHICH FOLLOW, THE ARGUMENT TYPES
C  ARE INDICATED AS FOLLOWS:
C               (DBL) IS FOR DOUBLE PRECISION
C               (INT) IS FOR INTEGER
C
C******************************************************************* 
      SUBROUTINE TDFRDR(IUNIT,NTR,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C   
C  TDFRDR READS THE DATA FROM TDF INTO THREE LABELLED COMMON BLOCKS
C  IN MEMORY: TDFCON, TDFSVB, AND TDFSPE.  TDFRDR MUST BE CALLED BY
C  THE APPLICATION CODE BEFORE ANY OTHER TDF ROUTINE CAN BE USED.
C  IT SHOULD BE CALLED ONLY ONCE.
C
C  NOTE: TDF MUST ALREADY HAVE BEEN OPENED BY THE APPLICATION
C  CODE AND ASSIGNED THE UNIT NUMBER IUNIT WHEN THIS CALL IS INVOKED. 
C
C  INPUT ARGUMENT
C    IUNIT (INT): LOGICAL UNIT NUMBER ASSIGNED TO TDF BY THE
C                   APPLICATION CODE.
C
C  OUTPUT ARGUMENTS
C      NTR (INT): NUMBER OF THERMONUCLEAR REACTIONS READ FROM TDF.
C     IERR (INT): FLAG INDICATING IF ANY ERROR OCCURRED IN THE CALL
C                   TO TDFRDR. A RETURN OF ZERO MEANS THAT TDF
C                   WAS READ AND THE MEMORY ALLOCATED CORRECTLY.
C                   A NON-ZERO RETURNED VALUE MEANS THAT TDF
C                   WAS NOT READ PROPERLY, AND IS FATAL.
C
C  VERSION DP3E,  12 NOVEMBER 1991
C
C**********************************************************************
C-------------STORAGE PARAMETERS--------------------------------
      PARAMETER (NM=32,NP=16,NR=16,NS=16)
      PARAMETER (NSV=NR*NM,NSP=NR*NS)
C  NR = MAX NO OF REACTIONS
C  NM = MAX NO OF MAXW-AVG POINTS PER REACTION
C  NS = MAX NO OF SPECTRA PER REACTION
C  NP = MAX NO OF POINTS IN EACH SPECTRUM
C  REFERENCE STORAGE INFORMATION:
C    MAX ALLOCATION OF INTEGERS IN COMMON BLOCK STORAGE
C      = 9+NR*(7+2*NS) = 633
C    DITTO FOR FLOATING POINT VARIABLES
C      = 2+NR*(7+3*NM+6*NS*(1+NP)) = 27762
C    MAX STORAGE (ITEM COUNT) FOR ALL COMMON BLOCKS
C      = 11+NR*(14+3*NM+8*NS+6*NS*NP) = 28395
C--------------TDFCON------------------------------------
C  DATA READ COUNT FOR EVERYTHING
      COMMON /TDFCON/ NUMR,NUMSVB,NUMSPE,NUMSP3,NUMSP4,
C  ACTUAL AND ALLOWED STORAGE, TEMP AND SVB MULT EXPONENTS
     X ISTOR,IDIM,ITHMUL,ISVMUL,
C  ITEM COUNTS FOR EACH REACTION
     X NSPECT(NR),NSIGV(NR),NPART(NR),
C  INTEGER IDENTIFIER FOR EACH REACTION IN THE ORDER
C  THEY APPEARED IN TDF.
     X IDREAC(NR),
C  DATA POINTERS FOR EACH REACTION (IPTR0 WORKS OFF THE
C  STORAGE WITH INDEX NSV AND IPTR3,4. DITTO FOR NSP)
     X IPTR0(NR),IPTR3(NR),IPTR4(NR),
C  LONG STORAGE ARRAYS (NUMBER OF POINTS IN SPECTRUM)
     X NUMPT3(NSP),NUMPT4(NSP)
      SAVE /TDFCON/
C-------------TDFSVB-------------------------------------
C  TEMP AND SIGMA-V BAR MULTS
      COMMON /TDFSVB/ THMUL,SVMUL,
C  MASS CONSTANTS FOR EACH REACTION
     X AM1(NR),AM2(NR),AM3(NR),AM4(NR),QVAL(NR),
C  LONG STORAGE ARRAYS (MAXWELL-AVERAGE STUFF, ALL REACTIONS)
     X TMPSVL(NSV),SIGVBL(NSV),EXTRL(NSV)
      SAVE /TDFSVB/
C-------------TDFSPE-------------------------------------
C  ZERO TEMP SPECTRUM ENERGIES, EACH REACTION
      COMMON /TDFSPE/ EN03(NR),EN04(NR),
C  LONG STORAGE ARRAYS FOR SPECTRA, ALL REACTIONS CONCATENATED
C  TEMP = TEMPERATURE
C  EN = ENERGY
C  SP = SPECTRUM PROFILE
C  PF = SPECTRUM PROFILE INTEGRAL
     X TEMP3(NSP),T3SQ(NSP),T3LN(NSP),
     X TEMP4(NSP),T4SQ(NSP),T4LN(NSP),
C  NOTE NR*NS*NP IS MAX STORAGE SET UP FOR EACH OF NEXT 6 ARRAYS
     X PF3(NSP,NP),EN3(NSP,NP),SP3(NSP,NP),
     X PF4(NSP,NP),EN4(NSP,NP),SP4(NSP,NP)
      SAVE /TDFSPE/
C--------------------------------------------------------
      CHARACTER*16 KEYWRD,LOKWRD
      CHARACTER*3 IHEAD
      CHARACTER*1 IH
      EQUIVALENCE (IH,IHEAD)
      DATA LOKWRD /'8A3E911112'/
C
C  INITIALIZE FLAGS AND COUNTERS
C
      DATA IERR1,IERR2,IERR3,IERR4 /1,1,1,1/
C  IR = REACTION COUNTER
C  K0 = MAXAV ITEM COUNTER
C  KI = PARTICLE I SPECTRUM COUNTER
C  JI = PARTICLE I SPECTRUM ITEM COUNTER
      DATA IR,K0,K3,K4,K5,J3,J4,J5 /0,0,0,0,0,0,0,0/
      NTR=0
      IERR=0
      ISTOR=0
      NUMR=0
      NUMSVB=0
      NUMSPE=0
      NUMSP3=0
      NUMSP4=0
      NSPNP=NSP*NP
C
C  SET UNITS FOR THE TEMPERATURE AND REACTION RATE ARGUMENTS
C  FOR CALLS TO SVBLU, MXAVLU AND SPECLU.  THE UNITS FOR THE
C  TEMPERATURES AND REACTION RATES IN TDF ARE MEV AND
C  BARN-CM/SEC RESPECTIVELY.  TO KEEP TO THESE UNITS FOR THE
C  SUBROUTINE CALL ARGUMENTS, SET ITHMUL = 0, ISVMUL = 0,
C  THMUL = 1.D0, AND SVMUL = 1.D0.  THE FOLLOWING ASSIGNMENTS
C  RESULT IN HAVING THE TEMPERATURE ARGUMENT IN KEV AND THE 
C  REACTION RATE ARGUMENT IN CM**3/SEC:
C
      ITHMUL=3
      ISVMUL=24
      THMUL=1.D-3
      SVMUL=1.D-24
C
C  LOOP TO READ TDF AND STUFF STORAGE
C
C  READ TDF LINE, CHECK FOR STATUS
  10  READ (IUNIT,6,IOSTAT=IOS) IHEAD
   6  FORMAT(A3)
      IF(IOS) 4,8,3
C  RETURN TO CALLING PROGRAM WITH FLAG FOR ANY ERRORS ENCOUNTERED
  12  IERR=IERR+16
      GO TO 4
   3  IERR=1
   4  IF((IERR1*IERR2*IERR3*IERR4).NE.1) IERR=IERR+2
      IF(NUMSVB.NE.K0.OR.NUMSP3.NE.J3.OR.NUMSP4.NE.J4
     X  .OR.NUMSPE.NE.K3.OR.NUMSPE.NE.K4) IERR=IERR+4
      IF(IR.EQ.0) IERR=IERR+8
C  REACTION COUNT
      NUMR=IR
      NTR=IR
C  MAX COUNT OF ITEMS STORED IN COMMON BLOCKS
      IDIM=11+14*NR+8*NSP+3*NSV+6*NSPNP
C  COUNT OF ITEMS ACTUALLY STORED
      ISTOR=11+14*IR+4*(K3+K4)+3*(K0+J3+J4)
      RETURN
C  LOOK FOR KEY SYMBOLS AND PROCESS ACCORDINGLY
   8  IF(IH.NE.'(') GO TO 10
C  START NEW REACTION DATA BLOCK
      IF(IHEAD.EQ.'(1)') THEN
C  CHECK THAT TDF VERSION IS THE ONE EXPECTED
C  IF NOT, SET FATAL ERROR FLAG AND RETURN
        READ (IUNIT,'(2X,A16)') KEYWRD
        IF(KEYWRD.NE.LOKWRD) THEN
          IERR=99
          RETURN
        ENDIF
C  TEST THAT THE PREVIOUS REACTION DATA
C  READ WAS PROPERLY COMPLETED
        IF((IERR1*IERR2*IERR3*IERR4).NE.1) GO TO 4
C  SET FLAGS AND POINTERS
        IERR1=0
        IERR2=0
        IERR3=0
        IERR4=0
        IR=IR+1
        IF(IR.GT.NR) GO TO 12
        IPTR3(IR)=K3
        IPTR4(IR)=K4
        IPTR0(IR)=K0
        GO TO 10
      ENDIF
C  PHYSICAL CONSTANTS, Q AND EX VALUE
      IF(IHEAD.EQ.'(2)') THEN
        READ (IUNIT,*) AMUMEV,FNSTNV,CELER,QVALUE,EXVAL
        QVAL(IR)=QVALUE
        IERR1=1
        GO TO 10
      ENDIF
C  REACTION MASSES
C ****  EXPAND LATER FOR 5TH MASS ****
      IF(IHEAD.EQ.'(3)') THEN
        READ (IUNIT,*) A1,A2,A3,A4,A5
        EN03(IR)=A4*QVAL(IR)/(A3+A4)
        EN04(IR)=A3*QVAL(IR)/(A3+A4)
        AM1(IR)=A1
        AM2(IR)=A2
        AM3(IR)=A3
        AM4(IR)=A4
        IERR2=1
        GO TO 10
      ENDIF
C  PARTICLE IDENTIFIERS, INTENDED ITEM COUNTS, REACTION LABEL NUMBER
      IF(IHEAD.EQ.'(4)') THEN
        READ (IUNIT,*) IAZ1,IAZ2,IAZ3,IAZ4,IAZ5,
     X  NSPECT(IR),NSIGV(IR),NPART(IR),IDREAC(IR)
        NUMSPE=NUMSPE+NSPECT(IR)
        NUMSVB=NUMSVB+NSIGV(IR)
        IERR3=1
        GO TO 10
      ENDIF
C  MAXWELL-AVERAGED DATA
      IF(IHEAD.EQ.'(5)') THEN
        DO 11 I=1,NSIGV(IR)
        READ (IUNIT,*,IOSTAT=IOS) A,B,C
        IF(A.EQ.1.0D+0) A = 1.000001D+0
        IF(IOS) 4,7,3
   7    K0=K0+1
        IF(K0.GT.NSV) GO TO 12
        TMPSVL(K0)=LOG(A)
        SIGVBL(K0)=LOG(B)
        EXTRL(K0)=LOG(C)
  11    CONTINUE
        IERR4=1
        GO TO 10
      ENDIF
C  SPECTRUM DATA HEADER
      IF(IHEAD.EQ.'(7)') THEN
        READ (IUNIT,*) THE,IPART,NPTS,AREA,SIGV
        IF(IPART.EQ.3) THEN
          K3=K3+1
          IF(K3.GT.NSP) GO TO 12
          TEMP3(K3)=THE
          T3SQ(K3)=SQRT(THE)
          T3LN(K3)=LOG(THE)
          NUMPT3(K3)=NPTS
          NUMSP3=NUMSP3+NPTS
          GO TO 10
        ENDIF
        IF(IPART.EQ.4) THEN
          K4=K4+1
          IF(K4.GT.NSP) GO TO 12
          TEMP4(K4)=THE
          T4SQ(K4)=SQRT(THE)
          T4LN(K4)=LOG(THE)
          NUMPT4(K4)=NPTS
          NUMSP4=NUMSP4+NPTS
          GO TO 10
        ENDIF
        GO TO 4
      ENDIF
C  SPECTRUM DATA
C  (ASSUMES THIS DIRECTLY FOLLOWS SPECTRUM HEADER)
      IF(IHEAD.EQ.'(8)') THEN
        DO 1 I=1,NPTS
        READ (IUNIT,*,IOSTAT=IOS) J,A,B,C
        IF(IOS) 4,9,3
   9    IF(J.NE.I) GO TO 3
        IF(C.LT.0.D0) C=1.D0+C
        IF(IPART.EQ.3) THEN
          J3=J3+1
          IF(J3.GT.NSPNP) GO TO 12
          EN3(K3,J)=A
          SP3(K3,J)=B
          PF3(K3,J)=C
        ENDIF
        IF(IPART.EQ.4) THEN
          J4=J4+1
          IF(J4.GT.NSPNP) GO TO 12
          EN4(K4,J)=A
          SP4(K4,J)=B
          PF4(K4,J)=C
        ENDIF
   1    CONTINUE
      ENDIF
      GO TO 10
      END
C*******************************************************************
      SUBROUTINE TDMASS(IREAC,BM1,BM2,BM3,BM4,BM5,Q,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  TDMASS RETURNS THE NUCLEAR (NOT THE ATOMIC) MASSES AND THE
C  Q-VALUE FOR A SPECIFIED REACTION IN TDF. THE MASS UNITS ARE
C  IN AMU WHERE THE ATOMIC MASS OF 12C IS EXACTLY 12 AMU.
C
C  INPUT ARGUMENT:
C     IREAC (INT): IDENTIFIER OF THE REACTION FOR WHICH
C                    MASSES AND Q-VALUE ARE DESIRED
C
C  OUTPUT ARGUMENTS:
C       BM1 (DBL): MASS OF LIGHTER REACTING PARTICLE (AMU)
C       BM2 (DBL): MASS OF HEAVIER REACTING PARTICLE (AMU)
C       BM3 (DBL): MASS OF LIGHTEST PRODUCT PARTICLE (AMU)
C       BM4 (DBL): MASS OF NEXT LIGHTEST PRODUCT PARTICLE (AMU).
C       BM5 (DBL): ZERO FOR TWO-BODY FINAL STATE REACTIONS, ELSE
C                   THE MASS OF THE REMAINING PRODUCT PARTICLE
C                   IN THREE-BODY FINAL STATE REACTIONS (AMU).
C         Q (DBL): Q-VALUE OF THE REACTION (MEV)
C      IERR (INT): FLAG INDICATING IF ANY ERROR OCCURRED
C                   DURING THE CALL TO TDMASS.  ZERO MEANS
C                   CALL COMPLETED CORRECTLY.
C
C**********************************************************************
      PARAMETER (NM=32,NP=16,NR=16,NS=16)
      PARAMETER (NSV=NR*NM,NSP=NR*NS)
      COMMON /TDFCON/ NUMR,NUMSVB,NUMSPE,NUMSP3,NUMSP4,
     X ISTOR,IDIM,ITHMUL,ISVMUL,
     X NSPECT(NR),NSIGV(NR),NPART(NR),
     X IDREAC(NR),
     X IPTR0(NR),IPTR3(NR),IPTR4(NR),
     X NUMPT3(NSP),NUMPT4(NSP)
      SAVE /TDFCON/
      COMMON /TDFSVB/ THMUL,SVMUL,
     X AM1(NR),AM2(NR),AM3(NR),AM4(NR),QVAL(NR),
     X TMPSVL(NSV),SIGVBL(NSV),EXTRL(NSV)
      SAVE /TDFSVB/
C
C  LOOK UP MASSES AND Q VALUE FOR REACTION LABEL IREAC
C
C  CLEAR OUTPUT AND ERROR FLAG
      BM1=0.D0
      BM2=0.D0
      BM3=0.D0
      BM4=0.D0
      BM5=0.D0
      Q=0.D0
      IERR=0
C  FIND REACTION INDEX IR
      DO 3 IR=1,NUMR
      IF(IREAC.EQ.IDREAC(IR)) GO TO 4
   3  CONTINUE
C  ERROR RETURN--NO REACTION FOUND
      IERR=502
      RETURN
C  REACTION FOUND, RETURN STUFF
   4  Q=QVAL(IR)
      BM1=AM1(IR)
      BM2=AM2(IR)
      BM3=AM3(IR)
      BM4=AM4(IR)
C  TEMPORARY EXPEDIENT TO CALCULATE FIFTH MASS IN A THREE-BODY
C  FINAL STATE REACTION
      BM5=BM1+BM2-BM3-BM4-Q/931.5012D0
      IF(BM5.LT..1D0) BM5=0.D0
      RETURN
      END      
C*******************************************************************
      SUBROUTINE MXAVLU(IREAC,TP,SV,E1,E2,E3,E4,E5,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  MXAVLU OBTAINS MAXWELLIAN-AVERAGED REACTION RATES AND KINETIC
C  ENERGIES OF THE REACTING AND FINAL PARTICLES FOR SPECIFIED
C  REACTION AND PLASMA TEMPERATURE.
C
C  INPUT ARGUMENTS
C     IREAC (INT): THE IDENTIFIER OF THE REACTION FOR WHICH
C                     REACTION RATES AND MAXWELLIAN-AVERAGED
C                     ENERGIES ARE DESIRED
C        TP (DBL): PLASMA TEMPERATURE (KEV)
C
C  OUTPUT ARGUMENTS
C        SV (DBL): MAXWELL-AVERAGED REACTION RATE (CM**3/SEC)
C        E1 (DBL): MAXWELLIAN-AVERAGED KINETIC ENERGY OF THE FIRST
C                    REACTING PARTICLE (MEV)
C        E2 (DBL): MAXWELLIAN-AVERAGED KINETIC ENERGY OF THE SECOND
C                    REACTING PARTICLE (MEV)
C        E3 (DBL): MAXWELLIAN-AVERAGED KINETIC ENERGY OF THE FIRST
C                    FINAL PRODUCT (MEV)
C        E4 (DBL): MAXWELLIAN-AVERAGED KINETIC ENERGY OF THE SECOND
C                    FINAL PRODUCT (MEV)
C        E5 (DBL): ZERO FOR 2-BODY FINAL STATE REACTION, ELSE 
C                    MAXWELLIAN-AVERAGED KINETIC ENERGY OF THE THIRD
C                    FINAL PRODUCT (MEV)
C      IERR (INT): FLAG INDICATING IF ANY ERROR OCCURRED
C                    DURING THE CALL TO MXAVLU.  ZERO MEANS
C                    CALL COMPLETED CORRECTLY
C
C  HIERARCHY:
C  MXAVLU CALLS INTERPOLATING ROUTINES TDCLAG AND TDCHLU
C
C  DESIGN:
C  MXAVLU IDENTIFIES THE DATA ARRAY INDEX FOR THE SPECIFIED REACTION,
C  AND DETERMINES THE FOUR SEQUENTIAL DATA TEMPERATURES WHICH
C  BRACKET THE USER'S REQUESTED TEMPERATURE.  IT CALLS TDCLAG TO 
C  INTERPOLATE THE REACTION RATE AT THE USER'S REQUESTED TEMPERATURE.
C  IT CALLS TDCHLU TO INTERPOLATE THE MEAN RELATIVE KINETIC ENERGY,
C  ER, OF THE TWO REACTING PARTICLES.  FROM THIS AND THE PLASMA
C  TEMPERATURE IT COMPUTES THE MAXWELLIAN-AVERAGED ENERGIES OF ALL
C  PARTICLES OF THE REACTION, NAMELY, E1 THROUGH E5.
C
C**********************************************************************
      PARAMETER (NM=32,NP=16,NR=16,NS=16)
      PARAMETER (NSV=NR*NM,NSP=NR*NS)
      COMMON /TDFCON/ NUMR,NUMSVB,NUMSPE,NUMSP3,NUMSP4,
     X ISTOR,IDIM,ITHMUL,ISVMUL,
     X NSPECT(NR),NSIGV(NR),NPART(NR),
     X IDREAC(NR),
     X IPTR0(NR),IPTR3(NR),IPTR4(NR),
     X NUMPT3(NSP),NUMPT4(NSP)
      SAVE /TDFCON/
      COMMON /TDFSVB/ THMUL,SVMUL,
     X AM1(NR),AM2(NR),AM3(NR),AM4(NR),QVAL(NR),
     X TMPSVL(NSV),SIGVBL(NSV),EXTRL(NSV)
      SAVE /TDFSVB/
C  CLEAR OUTPUT AND ERROR FLAG
      SV=0.D0
      E1=0.D0
      E2=0.D0
      E3=0.D0
      E4=0.D0
      E5=0.D0
      IERR=0
C  INPUT
      TIN=THMUL*TP
C  ERROR CHECK
      IF(TIN.LE.0.D0) THEN
        IERR=101
        RETURN
      ENDIF
C  FIND REACTION INDEX IR
      DO 3 IR=1,NUMR
      IF(IREAC.EQ.IDREAC(IR)) GO TO 4
   3  CONTINUE
C  ERROR RETURN--NO REACTION FOUND
      IERR=102
      RETURN
C  DATA LOOP LIMITS AND REACTION CONSTANTS
   4  ISTART=IPTR0(IR)+3
      IEND=IPTR0(IR)+NSIGV(IR)
      A1=AM1(IR)
      A2=AM2(IR)
      A3=AM3(IR)
      A4=AM4(IR)
      QVALUE=QVAL(IR)
C  LOG INPUT
      XIN=LOG(TIN)
C  ERROR RETURN IF INPUT BELOW TABLE
      IF(XIN.LT.(TMPSVL(ISTART-2)-1.D-3)) THEN
        IERR=101
        RETURN
      ENDIF
C  FIND AND SET UP LOG TEMP QUARTET
      DO 1 I=ISTART,IEND
      IF(XIN.LE.TMPSVL(I)) GO TO 2
   1  CONTINUE
C  ERROR RETURN--INPUT TEMP TOO HIGH
      IERR=103
      RETURN
   2  J=I
      IF(I.EQ.IEND) J=IEND-1
      X1=TMPSVL(J-2)
      X2=TMPSVL(J-1)
      X3=TMPSVL(J)
      X4=TMPSVL(J+1)
C  SIG V BAR LOOKUP
      F1=SIGVBL(J-2)
      F2=SIGVBL(J-1)
      F3=SIGVBL(J)
      F4=SIGVBL(J+1)
      CALL TDCLAG(X1,X2,X3,X4,F1,F2,F3,F4,XIN,YOUT,IERR)
      IF(IERR.NE.0) THEN
        IERR=104
        RETURN
      ENDIF
      SV=SVMUL*EXP(YOUT)
C  ER LOOKUP
      F1=EXTRL(J-2)
      F2=EXTRL(J-1)
      F3=EXTRL(J)
      F4=EXTRL(J+1)
      CALL TDCHLU(X1,X2,X3,X4,F1,F2,F3,F4,XIN,YOUT,IERR)
      IF(IERR.NE.0) THEN
        IERR=105
        RETURN
      ENDIF
      ER=EXP(YOUT)
C  CALCULATE OTHER ITEMS
      EC=1.5D0*TIN
      E1=(A2*ER+A1*EC)/(A1+A2)
      E2=(A1*ER+A2*EC)/(A1+A2)
      E3=(A4*(ER+QVALUE)+A3*EC)/(A3+A4)
      E4=(A3*(ER+QVALUE)+A4*EC)/(A3+A4)
C  TEMPORARY EXPEDIENT TO CALCULATE OUTPUT PARTICLE AVERAGED
C  ENERGIES IN THE THREE-PARTICLE FINAL STATE REACTION
C  T(T,2N)ALPHA (REACTION IDENTIFIER 3)
      IF(IREAC.EQ.3) THEN
        ETOT=E1+E2+11.332D0
        E3=ETOT*0.4731D0
        E4=ETOT*0.3640D0
        E5=ETOT*0.1629D0
      ENDIF
      RETURN
      END
C*******************************************************************
      FUNCTION SVBLU(IREAC,TP,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  SVBLU LOOKS UP THE MAXWELL-AVERAGED REACTION RATE FOR THE
C  SPECIFIED REACTION AND PLASMA TEMPERATURE.
C
C  INPUT ARGUMENTS
C     IREAC (INT): THE IDENTIFIER OF THE REACTION FOR WHICH
C                    MAXWELL-AVERAGED REACTION RATES ARE DESIRED.
C        TP (DBL): PLASMA TEMPERATURE (KEV)
C
C  OUTPUT ARGUMENT
C      IERR (INT): FLAG INDICATING IF ANY ERROR OCCURRED
C                    DURING THE CALL TO SVBLU.  ZERO MEANS
C                    CALL COMPLETED CORRECTLY.
C
C  RETURNED QUANTITY
C     SVBLU (DBL): MAXWELL-AVERAGED REACTION RATE (CM**3/SEC)
C
C  HIERARCHY:
C  SVBLU CALLS INTERPOLATING ROUTINE TDCLAG
C
C  DESIGN:
C  SVBLU IDENTIFIES THE DATA ARRAY INDEX FOR THE SPECIFIED REACTION,
C  AND DETERMINES THE FOUR SEQUENTIAL DATA TEMPERATURES WHICH
C  BRACKET THE USER'S REQUESTED TEMPERATURE.  IT CALLS TDCLAG TO 
C  INTERPOLATE THE REACTION RATE AT THE USER'S REQUESTED TEMPERATURE.
C
C**********************************************************************
      PARAMETER (NM=32,NP=16,NR=16,NS=16)
      PARAMETER (NSV=NR*NM,NSP=NR*NS)
      COMMON /TDFCON/ NUMR,NUMSVB,NUMSPE,NUMSP3,NUMSP4,
     X ISTOR,IDIM,ITHMUL,ISVMUL,
     X NSPECT(NR),NSIGV(NR),NPART(NR),
     X IDREAC(NR),
     X IPTR0(NR),IPTR3(NR),IPTR4(NR),
     X NUMPT3(NSP),NUMPT4(NSP)
      SAVE /TDFCON/
      COMMON /TDFSVB/ THMUL,SVMUL,
     X AM1(NR),AM2(NR),AM3(NR),AM4(NR),QVAL(NR),
     X TMPSVL(NSV),SIGVBL(NSV),EXTRL(NSV)
      SAVE /TDFSVB/
C  CLEAR OUTPUT AND ERROR FLAG
      IERR=0
      SVBLU=0.D0
C  INPUT
      TIN=THMUL*TP
C  ERROR CHECK     
      IF(TIN.LE.0.D0) THEN
        IERR=201
        RETURN
      ENDIF
C  FIND REACTION INDEX IR
      DO 3 IR=1,NUMR
      IF(IREAC.EQ.IDREAC(IR)) GO TO 4
   3  CONTINUE
C  ERROR RETURN--NO REACTION FOUND
      IERR=202
      RETURN
C  DATA LOOP LIMITS
   4  ISTART=IPTR0(IR)+3
      IEND=IPTR0(IR)+NSIGV(IR)
C  INPUT
      XIN=LOG(TIN)
C  ERROR RETURN IF INPUT BELOW TABLE
      IF(XIN.LT.(TMPSVL(ISTART-2)-1.D-3)) THEN
        IERR=201
        RETURN
      ENDIF
C  FIND TEMP INDEX AND SET UP LOG TEMP QUARTET
      DO 1 I=ISTART,IEND
      IF(XIN.LE.TMPSVL(I)) GO TO 2
   1  CONTINUE
C  ERROR RETURN--INPUT TEMP TOO HIGH
      IERR=203
      RETURN
   2  J=I
      IF(I.EQ.IEND) J=IEND-1
      X1=TMPSVL(J-2)
      X2=TMPSVL(J-1)
      X3=TMPSVL(J)
      X4=TMPSVL(J+1)
C  SIG V BAR LOOKUP
      F1=SIGVBL(J-2)
      F2=SIGVBL(J-1)
      F3=SIGVBL(J)
      F4=SIGVBL(J+1)
      CALL TDCLAG(X1,X2,X3,X4,F1,F2,F3,F4,XIN,YOUT,IERR)
      IF(IERR.NE.0) THEN
C  ERROR RETURN--OUT OF BOUNDS INPUT
        IERR=204
      ELSE
        SVBLU=SVMUL*EXP(YOUT)
      ENDIF
      RETURN
      END
C*******************************************************************
      SUBROUTINE SPECLU(IREAC,TP,IPART,PIN,EOUT,SOUT,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  SPECLU IS AN INVERSE PARTICLE SPECTRUM LOOKUP FOR THE
C  THERMONUCLEAR REACTIONS IN TDF.  GIVEN THE PLASMA TEMPERATURE
C  AND THE CUMULATIVE PROBABILITY OF EMISSION OF A GIVEN FINAL
C  STATE PARTICLE, IT RETURNS THE PARTICLE'S ENERGY AND NORMALIZED
C  SPECTRAL AMPLITUDE (I.E., RELATIVE PROBABILITY) CORRESPONDING
C  TO THIS BIVARIATE INPUT.
C
C  WARNING:  THIS ROUTINE IS CURRENTLY IMPLEMENTED FOR TWO-BODY
C  FINAL STATE REACTIONS ONLY.  IT WILL RETURN A NON-ZERO ERROR
C  FLAG AND ZERO ARGUMENTS IF ANY OTHER FINAL STATE REACTION
C  TYPE IS SELECTED.
C
C  INPUT ARGUMENTS
C       IREAC (INT): IDENTIFIER OF THE REACTION FOR WHICH
C                      SPECTRAL INFORMATION IS DESIRED
C          TP (DBL): PLASMA TEMPERATURE (KEV)
C       IPART (INT): IDENTIFIER OF PARTICLE FOR WHICH SPECTRAL
C                      INFORMATION IS DESIRED:
C                      3 FOR FIRST PRODUCT PARTICLE
C                      4 FOR SECOND PRODUCE PARTICLE
C                      5 FOR THIRD PRODUCT PARTICLE (NOT CURRENTLY
C                      IMPLEMENTED)
C         PIN (DBL): CUMULATIVE PROBABILITY (RANGE: 0.TO 1.,
C                      INCLUSIVE) FOR EMISSION OF PARTICLE
C                      SPECIFIED BY IPART AT GIVEN TEMPERATURE TP 
C
C  OUTPUT ARGUMENTS
C        EOUT (DBL): ENERGY (MEV) OF PARTICLE SPECIFIED BY IPART.
C        SOUT (DBL): NORMALIZED SPECTRUM AMPLITUDE (1./MEV)
C                      OF PARTICLE SPECIFIED BY IPART AND
C                      CORRESPONDING TO EOUT AND TP.  (TO GET THE
C                      ACTUAL SPECTRAL AMPLITUDE, MULTIPLY THIS
C                      QUANTITY BY THE CORRESPONDING REACTION RATE.)
C        IERR (INT): FLAG INDICATING IF ANY ERROR OCCURRED
C                      DURING THE CALL TO SPECLU.  ZERO MEANS
C                      CALL COMPLETED CORRECTLY
C
C  HIERARCHY:
C  SPECLU CALLS TWO INTERPOLATION ROUTINES TDCHLU AND TDMCHL AND
C  THE INVERSE LOOKUP ROUTINE TDSELU
C
C  DESIGN:
C  SPECLU FIRST CHECKS THE RANGE OF THE INPUT VALUE OF PIN
C  AND RETURNS SPECIAL VALUES OF EOUT AND SOUT IF PIN IS EXACTLY
C  ZERO OR ONE.  FOR INPUT VALUES OF PIN BETWEEN ZERO AND ONE,
C  SPECLU DETERMINES THE FOUR SEQUENTIAL DATA TEMPERATURES (AND
C  CORRESPONDING SPECTRAL DATA BLOCKS IN MEMORY) WHICH BRACKET
C  THE USER'S REQUESTED TEMPERATURE.  AT EACH OF THESE FOUR
C  TEMPERATURES IT CALLS TDSELU TO INTERPOLATE THE ENERGY (E)
C  AND SPECTRUM AMPLITUDE (S) AT THE USER'S SPECIFIED CUMULATIVE
C  PROBABILITY INPUT.  THE RESULTING FOUR VALUES EACH OF E AND
C  S ARE THEN INTERPOLATED AT THE USER'S SPECIFIED TEMPERATURE,
C  BY CALLS TO TDCHLU OR TDMCHL, AS APPROPRIATE, AND OUTPUT AS
C  EOUT AND SOUT, RESPECTIVELY
C
C**********************************************************************
      PARAMETER (NM=32,NP=16,NR=16,NS=16)
      PARAMETER (NSV=NR*NM,NSP=NR*NS)
      COMMON /TDFCON/ NUMR,NUMSVB,NUMSPE,NUMSP3,NUMSP4,
     X ISTOR,IDIM,ITHMUL,ISVMUL,
     X NSPECT(NR),NSIGV(NR),NPART(NR),
     X IDREAC(NR),
     X IPTR0(NR),IPTR3(NR),IPTR4(NR),
     X NUMPT3(NSP),NUMPT4(NSP)
      SAVE /TDFCON/
      COMMON /TDFSVB/ THMUL,SVMUL,
     X AM1(NR),AM2(NR),AM3(NR),AM4(NR),QVAL(NR),
     X TMPSVL(NSV),SIGVBL(NSV),EXTRL(NSV)
      SAVE /TDFSVB/
      COMMON /TDFSPE/ EN03(NR),EN04(NR),
     X TEMP3(NSP),T3SQ(NSP),T3LN(NSP),
     X TEMP4(NSP),T4SQ(NSP),T4LN(NSP),
     X PF3(NSP,NP),EN3(NSP,NP),SP3(NSP,NP),
     X PF4(NSP,NP),EN4(NSP,NP),SP4(NSP,NP)
      SAVE /TDFSPE/
      DIMENSION PFNODE(NP),ENNODE(NP),SPNODE(NP)
      DIMENSION XX(4),YY(4),GG(4),FF(4)
C  CLEAR OUTPUT AND ERROR FLAG
      EOUT=0.D0
      SOUT=0.D0
      IERR=0
C
C  CHECK AND PROCESS BIVARIATE INPUT
C
      TIN=THMUL*TP
      IF((PIN*(1.D0-PIN)).LT.0.D0) THEN
C  ERROR RETURN--INPUT PIN OUT OF BOUNDS
        IERR=300
        RETURN
      ENDIF
      IF(TIN.LE.0.D0) THEN
C  ERROR RETURN--INPUT TEMP NEG OR ZERO
        IERR=301
        RETURN
      ENDIF
C  FIND REACTION DATA BLOCK INDEX IR
      DO 6 IR=1,NUMR
      IF(IREAC.EQ.IDREAC(IR)) GO TO 7
   6  CONTINUE
C  ERROR RETURN--CAN'T FIND REACTION LABEL
      IERR=302
      RETURN
C  CHECK FINAL STATE PARTICLE COUNT
   7  IF(NPART(IR).NE.2) THEN
C  ERROR RETURN--REACTION DOES NOT HAVE TWO BODY FINAL STATE
C  (THIS WILL BE EXPANDED LATER TO ALLOW 3-BODY WHEN READY)
        IERR=303
        RETURN
      ENDIF
C  SET UP REACTION-RELATED INDICES AND CONSTANTS BY
C  CHOICE OF OUTPUT PARTICLE
      IF(IPART.EQ.3) THEN
        ETHR=.01D0*EN03(IR)
        ISTART=IPTR3(IR)+3
        IEND=IPTR3(IR)+NSPECT(IR)
        EINF=100.D0*EN03(IR)
        TBOTT=TEMP3(ISTART-2)
      ELSE
        IF(IPART.EQ.4) THEN
          ETHR=.01D0*EN04(IR)
          ISTART=IPTR4(IR)+3
          IEND=IPTR4(IR)+NSPECT(IR)
          EINF=100.D0*EN04(IR)
          TBOTT=TEMP4(ISTART-2)
        ELSE
C  ERROR RETURN--BAD OUTPUT PARTICLE INDEX
          IERR=304
          RETURN
        ENDIF
      ENDIF
C  CHECK IF INPUT TEMP BELOW TABLE
      IF(TIN.LT.(.999D0*TBOTT)) THEN
        IERR=301
        RETURN
      ENDIF
C
C  PSEUDO-LOOKUP FOR PIN = 0 OR 1
C
      IF(PIN.EQ.0.D0) THEN
        EOUT=0.D0
        SOUT=0.D0
        RETURN
      ENDIF
      IF(PIN.EQ.1.D0) THEN
        EOUT=EINF
        SOUT=0.D0
        RETURN
      ENDIF
C
C  BIVARIATE INVERSE SPECTRAL LOOKUP
C  DO LOOKUPS AT 4 BRACKET TEMPS, THEN
C  CUBIC HERMITE INTERP OVER THE 4 RESULTS
C
C  FIND 4 BRACKET TEMPS
      DO 1 I=ISTART,IEND
      IF(IPART.EQ.3) TT=TEMP3(I)
      IF(IPART.EQ.4) TT=TEMP4(I)
      IF(TIN.LE.TT) GO TO 5
  1   CONTINUE
C  ERROR RETURN--TEMP TOO HIGH
      IERR=305
      RETURN
  5   J=I
      IF(I.EQ.IEND) J=IEND-1
C  LOOK UP 4 E'S AND 4 S'S FOR 4 TEMPS
      N=0
      DO 2 K=J-2,J+1
      N=N+1
      IF(IPART.EQ.3) THEN
        DO 3 M=1,NUMPT3(K)
        PFNODE(M)=PF3(K,M)
        ENNODE(M)=EN3(K,M)
  3     SPNODE(M)=SP3(K,M)
        NNODE=NUMPT3(K)
        XX(N)=T3LN(K)
        YY(N)=T3SQ(K)
      ENDIF
      IF(IPART.EQ.4) THEN
        DO 4 M=1,NUMPT4(K)
        PFNODE(M)=PF4(K,M)
        ENNODE(M)=EN4(K,M)
  4     SPNODE(M)=SP4(K,M)
        NNODE=NUMPT4(K)
        XX(N)=T4LN(K)
        YY(N)=T4SQ(K)
      ENDIF
      CALL TDSELU(NNODE,PFNODE,ENNODE,SPNODE,PIN,E,S)
C  ERROR RETURN--NEGATIVE RESULTS
      IF(E.LE.0.D0.OR.S.LT.0.D0) THEN
        IERR=306
        RETURN
      ENDIF
      GG(N)=E
  2   FF(N)=LOG(S)
C  INTERPOLATE LOG S OVER LOG TEMPS
      X1=XX(1)
      X2=XX(2)
      X3=XX(3)
      X4=XX(4)
      F1=FF(1)
      F2=FF(2)
      F3=FF(3)
      F4=FF(4)
      XIN=LOG(TIN)
      CALL TDCHLU(X1,X2,X3,X4,F1,F2,F3,F4,XIN,FOUT,IERR)
      IF(IERR.NE.0) THEN
        IERR=307
        RETURN
      ENDIF
      SOUT=EXP(FOUT)
C  INTERPOLATE E OVER SQRT TEMPS
      YIN=SQRT(TIN)
      Y1=YY(1)
      Y2=YY(2)
      Y3=YY(3)
      Y4=YY(4)
      G1=GG(1)
      G2=GG(2)
      G3=GG(3)
      G4=GG(4)
      IF(G3.LT.ETHR) THEN
        IF(G2.LE.ETHR) THEN
          G1=LOG(G1)
          G2=LOG(G2)
          G3=LOG(G3)
          G4=LOG(G4)
          CALL TDCHLU(Y1,Y2,Y3,Y4,G1,G2,G3,G4,YIN,GOUT,IERR)
          IF(IERR.NE.0) THEN
            IERR=308
            RETURN
          ENDIF
          EOUT=EXP(GOUT)
        ELSE
          CALL TDMCHL(Y1,Y2,Y3,Y4,G1,G2,G3,G4,YIN,EOUT,IERR)
          IF(IERR.NE.0) THEN
            IERR=309
            RETURN
          ENDIF
        ENDIF
      ELSE
        CALL TDCHLU(Y1,Y2,Y3,Y4,G1,G2,G3,G4,YIN,EOUT,IERR)
        IF(IERR.NE.0) THEN
          IERR=310
          RETURN
        ENDIF
      ENDIF
      RETURN
      END
C*******************************************************************
      SUBROUTINE TDCLAG(X1,X2,X3,X4,F1,F2,F3,F4,X,FVAL,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  TDCLAG:  CUBIC LAGRANGE INTERPOLATION ROUTINE FOR MXAVLU AND SVBLU.
C  (INTERNAL ROUTINE NOT INTENDED TO BE CALLED DIRECTLY BY USER)
C
C**********************************************************************
      IERR=0
C  OUT OF BOUNDS TRAP
      IF(X.LT.X1.OR.X.GT.X4) THEN
        IERR=12
        RETURN
      ENDIF
C  LOOKUP
      FVAL=F1*(X-X2)*(X-X3)*(X-X4)/((X1-X2)*(X1-X3)*(X1-X4))+
     X     F2*(X-X1)*(X-X3)*(X-X4)/((X2-X1)*(X2-X3)*(X2-X4))+
     X     F3*(X-X1)*(X-X2)*(X-X4)/((X3-X1)*(X3-X2)*(X3-X4))+
     X     F4*(X-X1)*(X-X2)*(X-X3)/((X4-X1)*(X4-X2)*(X4-X3))
      RETURN
      END
C*******************************************************************
      SUBROUTINE TDSELU(NNODE,PFNODE,ENNODE,SPNODE,P,E,S)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  TDSELU:  SPECTRUM-ENERGY LOOKUP ROUTINE FOR USE BY TDFSPC;
C  INPUT P, RETURN E AND S ASSUMING P IS IN (0,1)
C  (INTERNAL ROUTINE NOT INTENDED TO BE CALLED DIRECTLY BY USER)
C
C**********************************************************************
      DIMENSION PFNODE(1),ENNODE(1),SPNODE(1)
C  LOWER EXTRAPOLATION INTERVAL
      IF(P.LE.PFNODE(1)) THEN
        X=P/PFNODE(1)
        A=PFNODE(1)/(ENNODE(1)*SPNODE(1))
        E=ENNODE(1)*X**A
        S=SPNODE(1)*X**(1.D0-A)
        RETURN
      ENDIF
C  UPPER EXTRAPOLATION INTERVAL
      IF(P.GE.PFNODE(NNODE)) THEN
        B=SPNODE(NNODE)/(1.D0-PFNODE(NNODE))
        S=B*(1.D0-P)
        E=ENNODE(NNODE)+LOG(SPNODE(NNODE)/S)/B
        RETURN
      ENDIF
C INTERPOLATION INTERVALS
      MID=(NNODE+1)/2
      ETHR=.5D0*ENNODE(MID)
      X=.5D0*LOG(P/(1.D0-P))
      DO 1 I=2,NNODE
      IF(P.LE.PFNODE(I)) GO TO 2
  1   CONTINUE
  2   ELEFT=ENNODE(I-1)
C  LEFT END
      P1=PFNODE(I-1)
      X1=.5D0*LOG(P1/(1.D0-P1))
      F1=ENNODE(I-1)
      S1=2.D0*P1*(1.D0-P1)/SPNODE(I-1)
C  RIGHT END
      P2=PFNODE(I)
      X2=.5D0*LOG(P2/(1.D0-P2))
      F2=ENNODE(I)
      S2=2.D0*P2*(1.D0-P2)/SPNODE(I)
C  SWITCH INTERP MODE IF CALLED FOR
      IF(ELEFT.LT.ETHR) THEN
        S1=S1/F1
        F1=LOG(F1)
        S2=S2/F2
        F2=LOG(F2)
      ENDIF
C  COEFS
      XD=1.D0/(X2-X1)
      SS=S1+S2
      SX=X2+X1
      SX2=X2**2+X1**2
      A=(SS-2.D0*(F2-F1)*XD)*XD**2
      B=.5D0*((S2-S1)*XD-3.D0*A*SX)
      C=.5D0*SS-1.5D0*A*SX2-B*SX
      D=.5D0*(F1+F2-A*(X2**3+X1**3)-B*SX2-C*SX)
C  OUTPUT VALUES
      E=D+X*(C+X*(B+X*A))
      ESL=C+X*(2.D0*B+X*3.D0*A)
      S=2.D0*P*(1.D0-P)/ESL
C  SWITCH INTERP MODE IF CALLED FOR
      IF(ELEFT.LT.ETHR) THEN
        E=EXP(E)
        S=S/E
      ENDIF
      RETURN
      END
C*******************************************************************
      SUBROUTINE TDCHLU(X1,X2,X3,X4,F1,F2,F3,F4,X,FVAL,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  TDCHLU: 4-POINT CUBIC HERMITE-LAGRANGE SLOPE-VALUE
C  INTERPOLATION ROUTINE FOR MXAVLU AND SPECLU.
C  (INTERNAL ROUTINE NOT INTENDED TO BE CALLED DIRECTLY BY USER)
C
C**********************************************************************
      IERR=0
C  OUT OF BOUNDS TRAP
      IF(X.LT.X1.OR.X.GT.X4) THEN
        IERR=10
        RETURN
      ENDIF
C  LOWER INTERVAL
      IF(X.LE.X2) THEN
      FVAL=F1*(X-X2)*(X-X3)/((X1-X2)*(X1-X3))+
     X     F2*(X-X1)*(X-X3)/((X2-X1)*(X2-X3))+
     X     F3*(X-X1)*(X-X2)/((X3-X1)*(X3-X2))
      RETURN
      ENDIF
C  UPPER INTERVAL
      IF(X.GE.X3) THEN
      FVAL=F2*(X-X3)*(X-X4)/((X2-X3)*(X2-X4))+
     X     F3*(X-X2)*(X-X4)/((X3-X2)*(X3-X4))+
     X     F4*(X-X2)*(X-X3)/((X4-X2)*(X4-X3))
      RETURN
      ENDIF
C  MIDDLE INTERVAL
      DX32=X3-X2
      XD32=1.D0/DX32
      DF43=(F4-F3)/(X4-X3)
      DF32=(F3-F2)*XD32
      DF21=(F2-F1)/(X2-X1)
      S2=(F3-F1+(2.D0*X2-X1-X3)*(DF32-DF21))/(X3-X1)
      S3=(F4-F2+(2.D0*X3-X2-X4)*(DF43-DF32))/(X4-X2)
      DX3=X-X3
      DX2=X-X2
      G2=F2*(3.D0*X2-2.D0*X-X3)-S2*DX2*DX32
      G3=F3*(3.D0*X3-2.D0*X-X2)+S3*DX3*DX32
      FVAL=(G3*DX2**2-G2*DX3**2)*XD32**3
      RETURN
      END
C*******************************************************************
      SUBROUTINE TDMCHL(X1,X2,X3,X4,F1,F2,F3,F4,X,FVAL,IERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**********************************************************************
C
C  TDMCHL:  MIXED 4-POINT CUBIC HERMITE-LAGRANGE SLOPE-VALUE
C  INTERPOLATION ROUTINE FOR SPECLU.  INTERPOLATION IS IN LINEAR SPACE
C  ON NODES 1, 2 AND 3  AND IN LOGARITHMIC SPACE ON NODES 2, 3 AND 4.
C  (INTERNAL ROUTINE NOT INTENDED TO BE CALLED DIRECTLY BY USER)
C
C**********************************************************************
      IERR=0
C  OUT OF BOUNDS TRAP
      IF(X.LT.X1.OR.X.GT.X4) THEN
        IERR=11
        RETURN
      ENDIF
C  LOWER INTERVAL
      IF(X.LE.X2) THEN
      FVAL=F1*(X-X2)*(X-X3)/((X1-X2)*(X1-X3))+
     X     F2*(X-X1)*(X-X3)/((X2-X1)*(X2-X3))+
     X     F3*(X-X1)*(X-X2)/((X3-X1)*(X3-X2))
      RETURN
      ENDIF
C  UPPER INTERVAL
      IF(X.GE.X3) THEN
      HVAL=LOG(F2)*(X-X3)*(X-X4)/((X2-X3)*(X2-X4))+
     X     LOG(F3)*(X-X2)*(X-X4)/((X3-X2)*(X3-X4))+
     X     LOG(F4)*(X-X2)*(X-X3)/((X4-X2)*(X4-X3))
      FVAL=EXP(HVAL)
      RETURN
      ENDIF
C  MIDDLE INTERVAL
      DX32=X3-X2
      XD32=1.D0/DX32
C  GET LINEAR SLOPE AT 2, CONVERT TO LOG SLOPE
      DF21=(F2-F1)/(X2-X1)
      DF32=(F3-F2)*XD32
      S2=(F3-F1+(2.D0*X2-X1-X3)*(DF32-DF21))/(X3-X1)
      S2=S2/F2
C  CONVERT F'S TO LOGS, WORK IN LOG F SPACE
      H2=LOG(F2)
      H3=LOG(F3)
      H4=LOG(F4)
      DF32=(H3-H2)*XD32
      DF43=(H4-H3)/(X4-X3)
      S3=(H4-H2+(2.D0*X3-X2-X4)*(DF43-DF32))/(X4-X2)
C  MIXED INTERP
      DX3=X-X3
      DX2=X-X2
      G2=H2*(3.D0*X2-2.D0*X-X3)-S2*DX2*DX32
      G3=H3*(3.D0*X3-2.D0*X-X2)+S3*DX3*DX32
      HVAL=(G3*DX2**2-G2*DX3**2)*XD32**3
      FVAL=EXP(HVAL)
      RETURN
      END
C*******************************************************************