      SUBROUTINE QRJ
C     ****      CALCULATE HEATING AND DISSOCIATION RATES
C     ****      CHAPMAN INTEGRAL FOR EACH SPECIES AT EACH GRID
C     ****      POINT IN NNO2, NNO, NNN2 (CALCULATED BY SUB CHAPMN)
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      include "crates.h"
      include "cmpbnd.h"
      parameter (lmax=59)
      COMMON/MASS/ RMN4S,RMN2D,RMNO,BRN2D,CEE
      DIMENSION EUVEFF(ZKMXP),SIGEUV(3,lmax),FEUV(lmax),
     1  RLMEUV(lmax),FSRC(15),SIGSRC(15),RLMSRC(15)
      DIMENSION SIGMAS(6,lmax)
      DIMENSION QUENCH(4)
      DIMENSION SFLUX(lmax),wave1(lmax),wave2(lmax)
      DIMENSION BRN2(lmax)
      DIMENSION BRO2(lmax)
C     ****
C     ****  TOTAL ABSORPTION COEFFICIENTS
C     ****
C     ****  SIGEUV(1,N) = SIGA(O2)
C     ****  SIGEUV(2,N) = SIGA(O)
C     ****  SIGEUV(3,N) = SIGA(N2)
C     ****
C     ****  TOTAL IONIZATION COEFFICIENTS
C     ****
C     ****  SIGMAS(1,N) = SIGI(O2)
C     ****  SIGMAS(2,N) = SIGI(O+(4S))
C     ****  SIGMAS(3,N) = SIGI(N2)
C     ****  SIGMAS(4,N) = SIGI(N)
C     ****  SIGMAS(5,N) = SIGI(O+(2D))
C     ****  SIGMAS(6,N) = SIGI(O+(2P))
C     ****
C     ****
C
      DATA RLMEUV/0.17250E-04, 0.16750E-04, 0.16250E-04, 0.15750E-04,
     +           0.15250E-04, 0.14750E-04, 0.14250E-04, 0.13750E-04,
     +           0.13250E-04, 0.12750E-04, 0.12250E-04, 0.12157E-04,
     +           0.11750E-04, 0.11250E-04, 0.10750E-04, 0.10250E-04,
     +           0.10319E-04, 0.10257E-04, 0.97500E-05, 0.97702E-05,
     +           0.92500E-05, 0.87500E-05, 0.82500E-05, 0.77500E-05,
     +           0.78936E-05, 0.77041E-05, 0.76515E-05, 0.72500E-05,
     +           0.70331E-05, 0.67500E-05, 0.62500E-05, 0.62973E-05,
     +           0.60976E-05, 0.57500E-05, 0.58433E-05, 0.55437E-05,
     +           0.52500E-05, 0.47500E-05, 0.46522E-05, 0.42500E-05,
     +           0.37500E-05, 0.36807E-05, 0.32500E-05, 0.30378E-05,
     +           0.30331E-05, 0.27500E-05, 0.28415E-05, 0.25630E-05,
     +           0.22500E-05, 0.17500E-05, 0.12500E-05, 0.75000E-06,
     +           0.41000E-06, 0.27500E-06, 0.19500E-06, 0.12000E-06,
     +           0.60000E-07, 0.30000E-07, 0.15000E-07/
      DATA(SIGEUV(1,N),N=1,lmax)/
     +             0.50, 1.50, 3.40, 6.00,10.00,13.00,
     +            15.00,12.00, 2.20, 0.30, 3.00, 0.01,
     +             0.30, 0.10, 1.00, 1.10, 1.00, 1.60,
     +            16.53, 4.00,15.54, 9.85,20.87,27.09,
     +            26.66,25.18,21.96,29.05,25.00,26.27,
     +            26.02,25.80,26.10,25.04,22.00,25.59,
     +            24.06,21.59,20.40,19.39,18.17,18.40,
     +            17.19,16.80,16.80,15.10,15.70,13.20,
     +            10.60, 7.10, 4.00, 1.18, 0.32, 0.10,
     +             1.02, 0.14, .024, .004, .0004/
      DATA(SIGEUV(2,N),N=1,lmax)/
     +              18 * 0.00,
     +             0.00, 0.00, 2.12, 4.18, 4.38, 4.23,
     +             4.28, 4.18, 4.18, 8.00,11.35,10.04,
     +            12.21,12.22,12.23,11.90,12.17,12.13,
     +            11.91,11.64,11.25,11.21, 9.64, 9.95,
     +             8.67, 7.70, 7.68, 6.61, 7.13, 6.05,
     +             5.30, 2.90, 1.60, 0.59, 0.16, 0.05,
     +             0.51, 0.07, .012, .002, .0002/
      DATA(SIGEUV(3,N),N=1,lmax)/
     +              18 * 0.00,
     +            36.16, 0.70,16.99,46.63,15.05,30.71,
     +            19.26,26.88,35.46,30.94,26.30,29.75,
     +            23.22,23.20,23.10,22.38,23.20,24.69,
     +            24.53,21.85,21.80,21.07,17.51,18.00,
     +            13.00,11.60,11.60,10.30,10.60, 9.70,
     +             8.00, 4.40, 1.90, 0.60, 0.24, 1.16,
     +             0.48, 0.09, .015, .003, .0003/
      DATA(SIGMAS(1,N),N=1,lmax)/
     +             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     +             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     +             0.00, 0.00, 0.00, 0.27, 0.00, 1.00,
     +            12.22, 2.50, 9.34, 4.69, 6.12, 9.39,
     +            11.05, 9.69, 8.59,23.81,23.00,22.05,
     +            25.94,25.80,26.10,25.04,22.00,25.59,
     +            24.06,21.59,20.40,19.39,18.17,18.40,
     +            17.19,16.80,16.80,15.10,15.70,13.20,
     +            10.60, 7.10, 4.00, 1.18, 0.32, 0.10,
     +             1.02, 0.14, .024, .004, .0004/
      DATA(SIGMAS(2,N),N=1,lmax)/
     +              20 * 0.00, 2.12, 4.18, 4.38, 4.23,
     +             4.28, 4.18, 4.18, 4.20, 4.91, 4.01,
     +             3.78, 3.79, 3.67, 3.45, 3.53, 3.52,
     +             3.45, 3.26, 3.15, 3.03, 2.51, 2.59,
     +             2.25, 1.93, 1.92, 1.65, 1.78, 1.51,
     +             1.38, 0.78, 0.46, 0.18, 0.05, 0.015,
     +             0.015, 0.02, 0.004, 0.0006, 0.00006/
      DATA(SIGMAS(3,N),N=1,lmax)/
     +              18 * 0.00,
     +             0.00, 0.00, 0.00, 0.00, 0.00,16.75,
     +            10.18,18.39,23.77,23.20,23.00,25.06,
     +            23.22,23.20,23.10,22.38,23.20,24.69,
     +            24.53,21.85,21.80,21.07,17.51,18.00,
     +            13.00,11.60,11.60,10.30,10.60, 9.70,
     +             8.00, 4.40, 1.90, 0.60, 0.24, 1.16,
     +             0.48, 0.09, .015, .003, .0003/
      DATA(SIGMAS(4,N),N=1,lmax)/
     +              22 * 0.00,
     +             1.00, 1.00, 1.00, 1.00, 1.10, 1.10,
     +             1.10, 1.20, 1.20, 1.20, 1.20, 1.20,
     +             1.10, 1.20, 1.15, 1.10, 1.00, 1.00,
     +             1.00, 0.70, 0.80, 0.65, 0.60, 0.60,
     +             0.50, 0.50, 0.40, 0.35, 0.25, 0.20,
     +             0.10, 0.10, 0.10, 0.05, 0.01,
     +             2*0.0/
      DATA(SIGMAS(5,N),N=1,lmax)/
     +              27 * 0.00, 3.80, 6.44, 5.52, 5.49,
     +             5.50, 5.50, 5.36, 5.48, 5.46, 5.36,
     +             5.24, 5.06, 4.71, 3.86, 3.98, 3.47,
     +             2.85, 2.84, 2.38, 2.64, 2.12, 1.86,
     +             0.99, 0.51, 0.19, 0.05, 0.015, 0.015,
     +             0.02, 0.004, 0.0006, 0.00006/
      DATA(SIGMAS(6,N),N=1,lmax)/
     +              29 * 0.00, 0.50, 2.93, 2.93, 3.06,
     +             3.09, 3.16, 3.15, 3.10, 3.14, 3.04,
     +             3.48, 3.28, 3.38, 2.95, 2.93, 2.92,
     +             2.58, 2.71, 2.42, 2.07, 1.13, 0.62,
     +             0.22, 0.06, 0.02, 0.02, 0.03, 0.004,
     +             0.0007, 0.00007/
      DATA(BRN2(N),N=1,lmax)/
     +              36 * 0.00,
     +             0.01, 0.04, 0.04, 0.03, 0.05, 0.05,
     +             0.15, 0.20, 0.20, 0.25, 0.32, 0.34,
     +             11 * 0.36/
      DATA(BRO2(N),N=1,lmax)/
     +              30 * 0.00,
     +             .025, .036, .045, .120, .155, .189,
     +             .230, 0.20, 0.20, 0.20, 0.23, 0.25,
     +             0.29, 16 * 0.33/
      DATA EUVEFF/29*0.05/
     2DO2/8.203E-12/,AD/5.435E-9/,BD/1.232E-6/,
     3ABAND/0.143/,BBAND/9.64E8/,CBAND/9.03E-19/
      DATA DO22/1.1407E-11/
      DATA QUENCH/7.E-11,5.E-11,3.1401E-12,9.1E-3/
      DATA ISTART/0/
      DATA E3/0.33/
      F107  = C(61)
      F107A = C(62)
      ISCALE = 0
C     ****
C     ****       THE FLAG ISCALE CHOOSES THE FLUX MODELING METHOD
C     ****       FROM THE SUBROUTINE SSFLUX. ISCALE=0 DOES SOLOMON'S,
C     ****       ISCALE=1 DOES A LIN. INTERP.,AND ISCALE=2 DOES TOBISKA'
C     ****
      IF(J.EQ.1) THEN
        CALL SSFLUX(ISCALE,F107,F107A,FLYA,SFLUX)
C       ****
C       ****     TRANSFER VALUES OF SFLUX TO APPROPRIATE SLOTS
C       ****
        DO 4 N = 16,lmax
          FEUV(N) = SFLUX(N)
    4   CONTINUE
        DO 5 N = 1,15
          FSRC(N)   = SFLUX(N)
    5   CONTINUE
      ENDIF
C     ****
C     ****      INITIALIZATION UPON FIRST ENTRY
C     ****
      IF(ISTART.EQ.0) THEN
        ISTART=1
        DO 1 M=1,3
          DO 1 N=1,lmax
          SIGEUV(M,N)=SIGEUV(M,N)*1.E-18
          SIGMAS(M,N)=SIGMAS(M,N)*1.E-18
    1   CONTINUE
        DO 2 N = 1,lmax
          SIGMAS(4,N) = SIGMAS(4,N)*1.E-17
          SIGMAS(5,N) = SIGMAS(5,N)*1.E-18
          SIGMAS(6,N) = SIGMAS(6,N)*1.E-18
    2   CONTINUE
        DO 3 N = 1,15
          RLMSRC(N) = RLMEUV(N)
          SIGSRC(N) = SIGEUV(1,N)
    3   CONTINUE
      ENDIF
C     ****
C     ****
C     ****     COMPUTE S14 = RSQ, S15 = RSP
C     ****
C     ****
C     ****     S1 = TAU(R),  S2 = TAU(Q)
C     ****
      NNO2K = NNO2
      NNOK = NNO
      NNN2K = NNN2
      DO 6 I = 1,LEN3
        S1(I,1) = SIGEUV(1,49)*F(I,NNO2K)+SIGEUV(2,49)*F(I,NNOK)+
     1    SIGEUV(3,49)*F(I,NNN2K)
        S2(I,1) = SIGEUV(1,20)*F(I,NNO2K)+SIGEUV(2,20)*F(I,NNOK)+
     1    SIGEUV(3,20)*F(I,NNN2K)
C       ****
C       ****     S3 = TAU(1), S4 = TAU(2), S5 = TAU(3)
C       ****
        S3(I,1) = 1.3*S1(I,1)
        S4(I,1) = 2.0*S1(I,1)
        S5(I,1) = 2.5*S1(I,1)
C       ****
C       ****     IF(TAU.GT.9.0) TAU = 9.0
C       ****
        S1(I,1) = CVMGP(S1(I,1),9.,9.-S1(I,1))
        S2(I,1) = CVMGP(S2(I,1),9.,9.-S2(I,1))
        S3(I,1) = CVMGP(S3(I,1),9.,9.-S3(I,1))
        S4(I,1) = CVMGP(S4(I,1),9.,9.-S4(I,1))
        S5(I,1) = CVMGP(S5(I,1),9.,9.-S5(I,1))
C       ****
C       ****     TAU(N) = EXP(-TAU(N)) FOR N = 1,3,1
C       ****
        S3(I,1) = EXP(-S3(I,1))
        S4(I,1) = EXP(-S4(I,1))
        S5(I,1) = EXP(-S5(I,1))
C       ****
C       ****     S6 = EXP(-TAU(R)), S7 = EXP(-TAU(Q))
C       ****
        S6(I,1) = EXP(-S1(I,1))
        S7(I,1) = EXP(-S2(I,1))
C       ****
C       ****     S14 = RSP, S15 = RSQ
C       ****
        S14(I,1) = S6(I,1)+2.*(S3(I,1)+S4(I,1)+S5(I,1))
        S15(I,1) = 1.5*S6(I,1)/(S14(I,1)+S2(I,1)/S1(I,1)*S7(I,1))
        S14(I,1) = 2.4*S6(I,1)/S14(I,1)
    6 CONTINUE
C     ****
C     ****
C     ****     S1 = PSI(O2), S2 = PSI(O), S3 = PSI(N2), S4 = PSI(N4S)
C     ****
C     ****
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NPN4SK=NJ+NPN4S
C     ****
C     ****      LEVELS 2 THRU KMAXP1
C     ****
      DO 7 I=1,LEN2
        S1(I,2)=.5*(F(I,NPS1K)+F(I,NPS1K+1))
        S2(I,2)=.5*(F(I,NPS2K)+F(I,NPS2K+1))
C       S4(I,2)=.5*(F(I,NPN4SK)+F(I,NPN4SK+1))
        S4(I,2)=0.
    7 CONTINUE
C     ****
C     ****      LEVEL 1
C     ****
      DO 8 I=1,LEN1
      S1(I,1) = .5*((B(I,1,1)+1.)*F(I,NPS1K)+B(I,1,2)*F(I,NPS2K)+
     1  FB(I,1))
      S2(I,1) = .5*(B(I,2,1)*F(I,NPS1K)+(B(I,2,2)+1.)*F(I,NPS2K)+
     1  FB(I,2))
C       S4(I,1) = F(I,NPN4SK)**1.5/SQRT(F(I,NPN4SK+1))
        S4(I,1) = 0.
    8 CONTINUE
C     ****
C     ****      NOW CALCULATE S3 = PSI(N2)
C     ****
      DO 9 I=1,LEN3
      S3(I,1)=1.-S1(I,1)-S2(I,1)
    9 CONTINUE
C     ****
C     ****
C     ****     CALCULATE IONIZATION CONTRIBUTIONS IN S8 - S13
C     ****
C     ****
C     ****       S8 = O2 IONIZATION
C     ****       S9 = O IONIZATION
C     ****       S10 = N2 IONIZATION
C     ****       S11 = N4S IONIZATION
C     ****       S12 = O2 IONIZATION * BRO2
C     ****       S13 = N2 IONIZATION * BRN2
C     ****
C     ****
C     ****     CONTRIBUTIONS TO Q, J AND QTEF FROM EUV
C     ****
C     ****      INITIALIZE Q, RJ AND QTEF
C     ****
      NQK = NQ
      NRJK = NRJ
      NQTEFK = NQTEF
      NQOP2PK = NQOP2P
      NQOP2DK = NQOP2D
      DO 10 I = 1,LEN3
        F(I,NQK) = 0.
        F(I,NRJK) = 0.
        F(I,NQTEFK) = 0.
        F(I,NQOP2PK) = 0.
        F(I,NQOP2DK) = 0.
   10 CONTINUE
C     ****
C     ****      INITIALIZE S8 - S13
C     ****
      DO 11 I=1,6*LEN3
        S13(I,1)=0.
   11 CONTINUE
C     ****
C     ****      SUMMATION OVER WAVE LENGTH
C     ****
      DO 12 N=16,lmax
C       ****
C       ****      S5=SUM(O2,O,N2)(SIGMA*CHAPMAN)
C       ****      S6=SUM(O2,O,N2)(SIGMA*PSI/RMASS)
C       ****      S7=SUM(O2,O,N2,N4S)(SIGMAS)
C       ****
C       ****      INITIALIZE S5, S6, S7
C       ****
        DO 13 I=1,3*LEN3
          S7(I,1)=0.
   13   CONTINUE
C       ****
C       ****      SUMMATION OVER O2, O, N2
C       ****
        DO 14 M=1,3
        K=(3-M)*KMAXP1+1
        NNK=NNO2+(M-1)*KMAXP1
        DO 14 I=1,LEN3
          S5(I,1)=S5(I,1)+SIGEUV(M,N)*F(I,NNK)
          S6(I,1)=S6(I,1)+SIGEUV(M,N)*S3(I,K)/RMASS(M)
          S7(I,1)=S7(I,1)+SIGMAS(M,N)
   14   CONTINUE
        DO 15 I = 1,LEN3
C         ****
C         ****     ADD CONTRIBUTIONS FROM N4S TO S7
C         ****
          S7(I,1) = S7(I,1)+SIGMAS(4,N)
C         ****
C         ****     S5 = F(LAMDA)*EXP(-TAU)
C         ****
          S5(I,1) = FEUV(N)*EXP(-S5(I,1))
C         ****
C         ****     ADD CONTRIBUTIONS TO Q AND RJE
C         ****
          F(I,NQK) = F(I,NQK)+C(60)/RLMEUV(N)*S5(I,1)*S6(I,1)
          F(I,NRJK) = F(I,NRJK)+S5(I,1)*S7(I,1)
          F(I,NQOP2PK) = F(I,NQOP2PK)+S5(I,1)*SIGMAS(6,N)*
     1                   S2(I,1)/RMASS(2)
          F(I,NQOP2DK) = F(I,NQOP2DK)+S5(I,1)*SIGMAS(5,N)*
     1                   S2(I,1)/RMASS(2)
   15   CONTINUE
C       ****
C       ****      ADD IONIZATION FOR O2, O, N2 TO S8, S9, S10
C       ****
        DO 16 M=1,3
          K=(3-M)*KMAXP1+1
          DO 16 I=1,LEN3
            S10(I,K)=S10(I,K)+SIGMAS(M,N)*S5(I,1)*S3(I,K)/RMASS(M)
   16   CONTINUE
        DO 17 I = 1,LEN3
C         ****
C         ****     ADD EFFECTIVE N2 IONIZATION TO QTEF
C         ****
          F(I,NQTEFK) = F(I,NQTEFK)+(SIGEUV(3,N)-SIGMAS(3,N))*S5(I,1)*
     1      S3(I,1)/RMASS(3)
C         ****
C         ****     ADD N4S IONIZATION TO S11
C         ****
          S11(I,1) = S11(I,1)+SIGMAS(4,N)*S5(I,1)*S4(I,1)/RMN4S
C         ****
C         ****     ADD IONIZATION CONTRIBUTIONS TO S12 AND S13
C         ****
          S12(I,1) = S12(I,1)+SIGMAS(1,N)*S5(I,1)*S1(I,1)/RMASS(1)*
     1      BRO2(N)
          S13(I,1) = S13(I,1)+SIGMAS(3,N)*S5(I,1)*S3(I,1)/RMASS(3)*
     1      BRN2(N)
   17   CONTINUE
   12 CONTINUE
C     ****
C     ****      MULTIPLY Q BY EFICIENCY FACTOR
C     ****
C     ****
      NQK=NQK-1
      DO 18 K=1,KMAXP1
        NQK=NQK+1
        DO 18 I=1,LEN1
          F(I,NQK)=F(I,NQK)*EUVEFF(K)*C(85)
   18 CONTINUE
C     ****
C     ****
C     ****     CALCULATE N*MBAR IN S5
C     ****
C     ****
C     ****       S6 = T(TOTAL,K)
C     ****
C     ****       LEVEL 1
C     ****
      NTK=NJ+NT+KMAX
      DO 19 I=1,LEN1
        S6(I,1)=F(I,NTK)+T0(1)
   19 CONTINUE
C     ****
C     ****      LEVELS 2 THRU KMAX
C     ****
      NTK=NJ+NT-1
      DO 20 K=2,KMAX
        NTK=NTK+1
        DO 20 I=1,LEN1
        S6(I,K)=.5*(F(I,NTK)+F(I,NTK+1))+T0(K)
   20 CONTINUE
C     ****
C     ****      LEVEL KMAXP1
C     ****
      NTK=NJ+NT+KMAX-1
      DO 21 I=1,LEN1
        S6(I,KMAXP1)=F(I,NTK)+T0(KMAXP1)
   21 CONTINUE
C     ****
C     ****      S7 = EXPS              (K)
C     ****
      DO 22 K=1,KMAX
        DO 22 I=1,LEN1
        S7(I,K)=C(87)*EXPS(K)
   22 CONTINUE
      DO 23 I=1,LEN1
        S7(I,KMAXP1)=C(86)*EXPS(KMAX)
   23 CONTINUE
C     ****
C     ****      S5 = N*MBAR            (K)
C     ****
      NMSK=NJ+NMS
      DO 24 I=1,LEN3
        S5(I,1)=C(81)*S7(I,1)*F(I,NMSK)/(C(84)*S6(I,1))
   24 CONTINUE
C     ****
C     ****
C     ****      CALCULATE CONTRIBUTIONS TO NQO2P, NQOP, NQN2P, NQNP
C     ****        AND NQTEF
C     ****
C     ****
      NQO2PK = NQO2P
      NQOPK = NQOP
      NQN2PK = NQN2P
      NQNPK = NQNP
      NQTEFK = NQTEF
      NQOP2PK = NQOP2P
      NQOP2DK = NQOP2D
      DO 25 I=1,LEN3
        F(I,NQO2PK) = F(I,NQO2PK)+(S8(I,1)*(1.+S15(I,1))-S12(I,1))*
     1    S5(I,1)
        F(I,NQN2PK) = F(I,NQN2PK)+(S10(I,1)*(1.+S14(I,1))-S13(I,1))*
     1              S5(I,1)
        F(I,NQNPK) = F(I,NQNPK)+(S11(I,1)+S13(I,1))*S5(I,1)
        F(I,NQTEFK) = F(I,NQTEFK)*S5(I,1)
        S4(I,1) = F(I,NQOP2PK)+F(I,NQOP2DK)+S9(I,1)+S12(I,1)
        F(I,NQOP2PK) = (F(I,NQOP2PK)+S4(I,1)*S14(I,1)*0.22)*S5(I,1)
        F(I,NQOP2DK) = (F(I,NQOP2DK)+S4(I,1)*S14(I,1)*0.24)*S5(I,1)
        F(I,NQOPK) = F(I,NQOPK)+(S9(I,1)+S12(I,1)+S4(I,1)*S14(I,1)
     1               *0.56)*S5(I,1)
   25 CONTINUE
C     ****
C     ****      CALCULATE NO IONIZATION AND ADD TO NQNOP
C     ****
      NQNOPK=NQNOP
      NNO2K=NNO2
      NPNOK=NJ+NPNO
C     ****      LEVEL 1
      DO 26 I=1,LEN1
        F(I,NQNOPK)=F(I,NQNOPK)+BETA9(I,1)*F(I,NPNOK)
     1    *S5(I,1)/RMNO
   26 CONTINUE
C     ****      LEVELS 2 THRU KMAXP1
      DO 27 I=LEN1+1,LEN3
        F(I,NQNOPK)=F(I,NQNOPK)+BETA9(I,1)*.5*(F(I,NPNOK)+F(I,NPNOK-1))
     1    *S5(I,1)/RMNO
   27 CONTINUE
      DO 28 K=1,KMAXP1
        FACTOR=C(85)*C(81)/C(57)*EXPS(1)*C(86)**(2*K-3)
        DO 28 I=1,LEN1
          S8(I,K)=FACTOR/((S1(I,K)/RMASS(1)+S2(I,K)/RMASS(2)+
     1      S3(I,K)/RMASS(3))*S6(I,K))
   28 CONTINUE
      DO 29 I=1,LEN3
        S8(I,1)=S8(I,1)*(QUENCH(1)*S3(I,1)/RMASS(3)+
     1    QUENCH(2)*S1(I,1)/RMASS(1))
        S8(I,1)=QUENCH(3)*S8(I,1)/(QUENCH(4)+S8(I,1))
   29 CONTINUE
      NQK=NQ
      NRJK=NRJ
C     ****
C     ****      S7=SUM OVER WAVE LENGTH(SIGMA*F*EXP(-SIGMA*CHAPMAN)*
C     ****      (HC/LAMDA-DO2))
C     ****      INITIALIZE S7
C     ****
      DO 30 I=1,LEN3
        S7(I,1)=0.
   30 CONTINUE
C     ****
C     ****      SUMMATION OVER WAVE LENGTH
C     ****
      NNO2K=NNO2
      DO 31 N=1,15
        DO 31 I=1,LEN3
        S9(I,1)=SIGSRC(N)*FSRC(N)*EXP(-SIGSRC(N)*F(I,NNO2K))
        S7(I,1)=S7(I,1)+S9(I,1)*
     1    (C(60)/RLMSRC(N)-DO22+S8(I,1))
C       ****
C       ****      UPDATE RJ
C       ****
        F(I,NRJK)=F(I,NRJK)+S9(I,1)
   31 CONTINUE
C     ****
C     ****      UPDATE Q
C     ****
      DO 32 I=1,LEN3
        F(I,NQK)=F(I,NQK)+S7(I,1)*C(85)*S1(I,1)/RMASS(1)
   32 CONTINUE
C     ****      CONTRIBUTIONS FROM SCHUMANN RUNGE BANDS
C     ****      S7=P3*F
      NNO2K=NNO2
      DO 33 I=1,LEN3
        S8(I,1)=SQRT(F(I,NNO2K))
        S7(I,1)=CVMGP(1./(ABAND*F(I,NNO2K)+BBAND*S8(I,1)),CBAND,
     1    F(I,NNO2K)-1.E18)*(1.+0.11*(C(61)-65.)/165.)
        F(I,NQK)=F(I,NQK)+S7(I,1)*C(85)*S1(I,1)/RMASS(1)*E3
        F(I,NRJK)=F(I,NRJK)+S7(I,1)/DO2
   33 CONTINUE
      RETURN
      END
C
