c
C IRIFLIP.for 
C
C 2012.00 10/05/11 IRI-2012: bottomside B0 B1 model (SHAMDB0D, SHAB1D),
C 2012.00 10/05/11    bottomside Ni model (iriflip.for), auroral foE
C 2012.00 10/05/11    storm model (storme_ap), Te with PF10.7 (elteik),
C 2012.00 10/05/11    oval kp model (auroral_boundary), IGRF-11(igrf.for), 
C 2012.00 10/05/11    NRLMSIS00 (cira.for), CGM coordinates, F10.7 daily
C 2012.00 10/05/11    81-day 365-day indices (apf107.dat), ap->kp (ckp),
C 2012.00 10/05/11    array size change jf(50) outf(20,1000), oarr(100).
C 2012.01 12/12/11 Deleted ALT_RATES (not used)
C 2012.02 01/04/12 Deleted FINDAP,READAP,CONV_DATE,GET_DATA,RATCHK (not used)
C 2012.02 01/04/12 Deleted BRACE,ACTUAL_DAY,EPHEM SOLDEC,TFILE,RUN_ERROR (not used)
C 2012.02 01/04/12 COP2D: 99 FOMRAT ',' missing; commented out all WRITEs
C 2012.03 07/17/14 COP4S: NPLUS=0; PR(13)=0.0 ------------------------- A Shabanloui
C 2016.01 09/08/16 Main: NEWTON replaced by iteration procedure ------- B Gustavsson
C 2016.02 05/07/18 Added array PRV11 to DATA statement ---------------- K. Knight
C 2020.01 03/05/24 SECIPRD: added IMAX=0 at beginning ----------------- I. Alexeyevich
C 2020.07 07/22/25 SCOLUM: IF(Y.GE.8) ERFY2=F/(G+Y) instead of GT ----- K. Johnston
C****************************************************************************************
C subroutines for IDC model
C
C includes: main subroutine CHEMION and the following subroutines and functions
C           KEMPPRN.FOR: CN2D, CNO, CN4S, CN2PLS, CNOP, CO2P, COP4S, COP2D, COP2P,
C                        CNPLS, CN2A, CN2P, CNOPV
C           RATES.FOR:   RATS 
C           PESIMP.FOR:  SECIPRD, FLXCAL, FACFLX, SIGEXS, TXSION, OXRAT, T_XS_N2, 
C                        T_XS_OX, OXSIGS
C           RSPRIM.FOR:  PRIMPR, SCOLUM, PARAMS, PROBS, PROBN2, YLDISS, PROBO2, 
C                        SCHUMN, FACEUV, FACSR,  
C  
C turn on printout of intermediate quantities with JPRINT=1 also in PARAMS, PROBS, 
C PROBN2, YLDISS, and PROBO2 with ISW=1.
C 
C Richards, P. G., D. Bilitza, and D. Voglozin (2010), Ion density calculator (IDC): 
C    A new efficient model of ionospheric ion densities, Radio Sci., 45, RS5007, 
C    doi:10.1029/2009RS004332.
C
C Fortran code written by Phil Richards, George Mason University, Fairfax, VA, USA
C
C****************************************************************************************
C
C
C:::::::::::::::::::::::::::: CHEMION :::::::::::::::::::::::::::
C..... This routine was written by Phil Richards April 2010. This version was
C..... modified in April 2011. 
C..... It takes the specified input electron density and returns O+, O2+, NO+,
C..... N2+, N+, NO, and N(2D) densities. These densities generally agree well
C..... with Atmosphere Explorer and FLIP model densities.
C..... In this version all the densities except O+ are calculated from chemical
C.....  equilibrium. The densities are normalized so that the total ion density 
C..... matches the input electron density.
C..... This version has two modes. If the variable USER_OPLUS is positive it is used
C..... to specify the O+ density. If USER_OPLUS is negative, O+ is calculated from
C..... chemical equilibrium. In both cases, all the ion densities are normalized
C..... to the input electron density (Ne). Thus the O+ density may not match
C..... exactly, USER_OPLUS.
C..... N+ generally agrees well with AE-C data and the FLIP model during the day
C..... up to ~500 km, but is inaccurate at night due to diffusion.
C..... The NO densities can either be user specified or calculated by the model.
C..... NO will be very good except below about 130 km where it will be 
C..... underestimated due to neglect of diffusion. There is an artificial 
C..... floor on the NO density to prevent it from getting too low below 130 km.
C..... H+ and He+ are only good during the daytime below ~450 km.
C..... The EUVAC model is used for solar EUV irradiances
      SUBROUTINE CHEMION(JPRINT,  !.. Input: Turn file output on or off
     >                      ALT,  !.. Input: Altitude(km)
     >               F107,F107A,  !.. Input: Solar activity indices
     >                 TE,TI,TN,  !.. Input: Electron, ion and neutral temperatures
     >       OXN,O2N,N2N,HEN,HN,  !.. Input: O, O2, N2, He, and H densities (cm-3)
     >                  USER_NO,  !.. Input: User specified NO density (cm-3)
     >                      N4S,  !.. Input: N4S should be 0.5*MSIS N density (cm-3)
     >                       NE,  !.. Input: electron density (cm-3)
     >               USER_OPLUS,  !.. Input: User specified O+ density (cm-3) -1.0=off
     >                     SZAD,  !.. Input: LT(hrs), UT(sec) and solar zenith angle(D)
     >     OXPLUS,O2PLUS,NOPLUS,  !.. OUTPUT: O+, O2+, NO+ densities (cm-3)
     >             N2PLUS,NPLUS,  !.. OUTPUT: N2+ and N+ densities (cm-3)
     >                  NNO,N2D,  !.. OUTPUT: NO and N(2D) density (cm-3)
     >                    ITERS)  !.. OUTPUT: # of iterations to converge
      IMPLICIT NONE
      INTEGER JPRINT        !.. Turns on printing of production and loss
      INTEGER I,J,K,ITERS   !.. loop control variables
      INTEGER IRATS         !.. Switch for different rates
      INTEGER ITS,JITER     !.. Variables for Newton procedure
      REAL TE,TN,TI         !.. Electron and ion temperatures
      !.. Geophysical parameters
      REAL F107,F107A,ALT,SZAD
      !.. Measured H+, He+, O+, N2+, NO+, O2+, N+, RPA ion density
      REAL HEPLUS,OXPLUS,N2PLUS,NOPLUS,O2PLUS,NPLUS,USER_Oplus
      !.. O2,O,N2,NO,He,N4S, user specified NO
      REAL O2N,OXN,N2N,NNO,HEN,N4S,HN,USER_NO
      !.. Ne, N(2P),N(2D),O+(2P),O+(2D) densities
      REAL NE,N2P,N2D,OP2D,OP2P
      !.. Total (photon & photoel) production rates O+(4S),O+(2P),O+(2D),O2+
      REAL TPROD1,PDISOP,TPROD2,TPROD3,TPROD5
      !.. Total Production rates from all sources for NO+, O2+, 
      REAL TPNOP,O2PPROD
      !.. Production rates hv(e*)+N2->N+, hv+N->N+, Lyman-a -> NO+ 
      REAL DISNP,PHOTN,PLYNOP
      REAL PSEC                     !.. generic PE production
      REAL RTS(99)                  !.. Reaction rates array
      REAL SECPN2PLUS,EUVN2PLUS     !.. N2+ total production
      REAL H,DEX,FEX(2)             !.. used in Newton solver
      REAL SUMIONS                  !.. Sum of the major ions
      REAL PNO,LNO,PDNOSR           !.. Production and loss of NO
      REAL N2A                      !.. N2(A) density    
      REAL DISN2D,UVDISN,PN2D,LN2D  !.. Production and loss of N(2D)
      REAL N2APRD                   !.. PE production rate of N2(A)
      REAL PN4S,LN4S,DISN4S         !.. Production and loss of N(4S)
      REAL COXPLUS                  !.. Chemical equilibrium O+
      !.. various ionization and excitation rates by EUV and PE
      REAL EUVION,PEXCIT,PEPION,OTHPR1,OTHPR2,PRHEP
      REAL SUMSAVE                  !.. saved sum of ions for convergence
      COMMON/EUVPRD/EUVION(3,12),PEXCIT(3,12),PEPION(3,12),OTHPR1(6)
     >   ,OTHPR2(6)

      !.. initialize parameters
      DATA K/0/
      DATA PNO,LNO,PDNOSR,PLYNOP,N2A/5*0.0/
      DATA DISN2D,UVDISN/0.0,0.0/

      JITER=0      !.. Counts the number of Newton iterations
      N2P=0.0      !.. N(2P) density, not calculated here

      CALL RATS(0,TE,TI,TN,RTS)  !.. Get the reaction rates

      !.. PRIMPR calculates solar EUV production rates. 
      CALL PRIMPR(1,ALT,OXN,N2N,O2N,HEN,SZAD*0.01745,TN,F107,F107A,N4S)

      !.. Calculate secondary Production from photoelectrons
      CALL SECIPRD(ALT,SZAD,F107,F107A,TE,TN,OXN,O2N,N2N,
     >   NE,N2APRD)

      UVDISN=OTHPR1(1)                   !.. EUV dissociation rate of N2
      DISNP= EUVION(3,4)+EUVION(3,5)+EUVION(3,6)
     >   +0.1*(PEPION(3,1)+PEPION(3,2)+PEPION(3,3))  !.. Rydberg diss       
     >   +PEPION(3,4)+PEPION(3,5)+PEPION(3,6)       
      DISN2D=2.0*PEPION(3,1)+OTHPR2(3)
      DISN4S=2.0*PEPION(3,1)+OTHPR2(3)
      PRHEP=OTHPR1(2)                          !.. He+ photoionization 

      !.. initialize variables to avoid using left over values
      HEPLUS=0.0
      OXPLUS=0.0
      N2PLUS=0.0
      NOPLUS=0.0
      O2PLUS=0.0
      N2P=0.0
      N2D=0.0
      OP2D=0.0
      OP2P=0.0
      N2A=0.0
      SUMSAVE=0.0

      K=K+1        !.. If K=1 print headers in files

      !.. These species don't need to be iterated because they are at 
      !.. the top of the food chain
      !.. O+(2P) Calculate and print densities, production, loss
      PSEC=PEPION(1,3)           !.. Photoelectron production
      TPROD3=EUVION(1,3)+PSEC    !.. Add EUV and photoelectrons
      CALL COP2P(JPRINT,7,K,ALT,RTS,OXN,O2N,N2N,NE,OP2P,TPROD3,PSEC
     >    ,HEPLUS,N4S,NNO,TE)

      !.. O+(2D) Calculate and print densities, production, loss
      PSEC=PEPION(1,2)           !.. Photoelectron production
      TPROD2=EUVION(1,2)         !.. EUV
      CALL COP2D(JPRINT,8,K,ALT,RTS,OXN,O2N,N2N,NE,OP2D,TPROD2,OP2P
     >    ,HEPLUS,N4S,NNO,PSEC)

      !.. O+(4S) Calculate and print densities, production, loss. 
      TPROD1=EUVION(1,1)
      PDISOP=EUVION(2,4)+EUVION(2,5)+PEPION(2,4)+PEPION(2,5)
      CALL COP4S(JPRINT,4,K,ALT,RTS,OXN,O2N,N2N,NE,COXPLUS,TPROD1,OP2D
     >    ,OP2P,PEPION(1,1),PDISOP,N2PLUS,N2D,NNO,1.0,HEPLUS)

      !.. Make sure chemical O+ is not greater than Ne
      IF(COXPLUS.GT.NE) COXPLUS=NE 

      !.. Choose either user specified or chemical O+
      IF(USER_OPLUS.GT.0)  THEN
        !.. This average smooths out bumps in the input O+
        OXPLUS=(USER_OPLUS+COXPLUS)/2
      ELSE
        !.. Alternative chemical equilibrium O+ calculation
        OXPLUS=COXPLUS  
      ENDIF

      !.. N2(A) is used in calculating NO density
      CALL CN2A(JPRINT,27,K,ALT,RTS,OXN,O2N,N2N,NE,N2A,N2APRD,0.0,
     >     0.0,0.0)

      !.. Iterate through chemistry to improve results
      DO ITERS=1,5
        !.. N2+ Calculate and print densities, production, loss. 
        CALL CN2PLS(JPRINT,9,K,ALT,RTS,OXN,O2N,N2N,NE,N2PLUS,EUVION(3,1)
     >   ,EUVION(3,2),EUVION(3,3),PEPION(3,1),PEPION(3,2),PEPION(3,3)
     >   ,OP2D,OP2P,HEPLUS,NPLUS,NNO,N4S)

        !.. N(2D) Calculate and print densities, production, loss. 
        CALL CN2D(JPRINT,16,K,ALT,RTS,OXN,O2N,N2N,NOPLUS,NE,PN2D,LN2D
     >    ,N2PLUS,DISN2D,UVDISN,NPLUS,N2P,N2D,OXPLUS,NNO,N2A)
        N2D=PN2D/LN2D

        !.. N+ Calculate and print densities, production, loss. 
        PHOTN=OTHPR2(3)  !.. N+ photo production
        CALL CNPLS(JPRINT,10,K,ALT,RTS,OXN,O2N,N2N,NE,DISNP,NPLUS,
     >    OXPLUS,N2D,OP2P,HEPLUS,PHOTN,O2PLUS,N4S,OP2D,N2PLUS,NNO)

        !.. NO Calculate and print densities, production, loss.
        CALL CNO(JPRINT,15,K,ALT,RTS,OXN,O2N,N2N,NE,PNO,LNO
     >    ,N2D,N4S,N2P,NNO,O2PLUS,OXPLUS,OTHPR2(2),OTHPR2(1),N2A,NPLUS)
        NNO=PNO/LNO     !.. NO chemical equilibrium density

        !.. Set a floor on NO density, which is needed below ~150 km at night 
        IF(NNO.LT.1.0E8*EXP((100-ALT)/20)) NNO=1.0E8*EXP((100-ALT)/20)
        IF(USER_NO.GT.1.0) NNO=USER_NO  !.. substitute user specified value
        IF(NNO.GT.1.5E8) NNO=1.5E8      !.. Don't let NO get too big

        !.. NO+ Calculate and print densities, production, loss. 
        CALL CNOP(JPRINT,11,K,ALT,RTS,OXN,O2N,N2N,NE,TPNOP,NOPLUS,OXPLUS
     >    ,N2PLUS,O2PLUS,N4S,NNO,NPLUS,N2P,PLYNOP,1.0,N2D,OP2D)

        !.. O2+ Calculate and print densities, production, loss. 
        !.. EUV + PE production
        TPROD5=EUVION(2,1)+EUVION(2,2)+EUVION(2,3)+PEPION(2,1)+
     >       PEPION(2,2)+PEPION(2,3)
        CALL CO2P(JPRINT,12,K,ALT,RTS,OXN,O2N,N2N,NE,O2PPROD
     >    ,O2PLUS,TPROD5,OXPLUS,OP2D,N2PLUS,NPLUS,N4S,NNO,OP2P)

        SUMIONS=OXPLUS+NOPLUS+O2PLUS+NPLUS+N2PLUS

        !.. Chemical equilibrium densities are normalized to the input NE 
        !.. and return.
        IF(ITERS.EQ.5.OR.ABS(SUMSAVE-SUMIONS)/SUMIONS.LT.0.01) THEN
          OXPLUS=OXPLUS*NE/SUMIONS
          NOPLUS=NOPLUS*NE/SUMIONS
          O2PLUS=O2PLUS*NE/SUMIONS
          N2PLUS=N2PLUS*NE/SUMIONS
          NPLUS=NPLUS*NE/SUMIONS
          RETURN
        ENDIF
        SUMSAVE=SUMIONS  
      ENDDO

      RETURN
      END
C
C
C:::::::::::::::::::::: KEMPRN.FOR ::::::::::::::::::::::::::::::::::::::::::::::::
C..... This file contains the chemistry routines for ions and neutrals
C..... First N(2D)
      SUBROUTINE CN2D(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NOP,NE,P1,L1
     > ,N2PLS,DISN2D,UVDISN,NPLUS,N2P,N2D,OPLS,NNO,N2A)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=NOP*NE*RTS(50)
      PR(2)=N2PLS*NE*RTS(32)*RTS(11)
      PR(3)=N2PLS*ON*RTS(10)
      PR(4)=DISN2D
      PR(5)=RTS(63)*UVDISN
      PR(6)=RTS(65)*NPLUS*O2N
      PR(7)=N2P*RTS(57)
      PR(8)=RTS(27)*N2A*ON
      LR(1)=ON*RTS(15)
      LR(2)=O2N*RTS(16)
      LR(3)=NE*RTS(8)
      LR(4)=OPLS*RTS(29)
      LR(5)=RTS(61)
      LR(6)=RTS(41)*NNO
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)
      L1=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)
C....... EF is used to convert production rates to volume emission rates
      EF=1.0
C....... This line in for volume emission rates
C...      EF=RTS(61)*0.76/L1
      IF(JPT.EQ.1.AND.JPR.GT.0.AND.INT(EF+0.1).NE.1) WRITE(I,189)
      IF(JPT.EQ.1.AND.JPR.GT.0.AND.INT(EF+0.1).EQ.1) WRITE(I,191)
 189  FORMAT(/2X,'N(2D)',25X,'EMISSION',28X,':',20X,'Loss rate')
 191  FORMAT(/2X,'N(2D)',25X,'Production',36X,':',20X,'Loss rate')
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,96)
 96   FORMAT(2X,'ALT   [N2D]   NO++e   N2++e   N2++O    e+N2   hv+N2'
     >  ,3X,'N++O2   N(2P)   N2A+O    +O     +O2      +e     +O+'
     >  ,5X,'RAD     +NO')
      IF(JPR.GT.0) WRITE(I,7) Z,P1/L1,(PR(K)*EF,K=1,8)
     > ,(LR(K)*N2D,K=1,6)
      RETURN
 7    FORMAT(F6.1,1P,22E8.1)
      END
C:::::::::::::::::::::::::::::: NO ::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE CNO(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1,L1
     > ,N2D,N4S,N2P,NNO,O2P,OPLS,PDNOSR,PLYNOP,N2A,NPLUS)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=RTS(16)*O2N*N2D
      PR(2)=RTS(7)*O2N*N4S
      PR(3)=RTS(38)*N2P*O2N
      PR(4)=RTS(27)*N2A*ON
      PR(5)=RTS(22)*NPLUS*O2N          !.. Fox
      LR(1)=RTS(9)*N4S
      LR(2)=RTS(23)*O2P
      LR(3)=RTS(24)*OPLS
      LR(4)=RTS(41)*N2D
      LR(5)=PDNOSR
      LR(6)=PLYNOP
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)
      L1=LR(1)+LR(2)+LR(3) + LR(4) + (LR(5)+LR(6))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,192)
 192  FORMAT(/2X,'NO',17X,'PRODUCTION',20X,':',10X,'LOSS RATES'/
     > ,4X,'ALT',3X,'[NO]',5X,'[NO]c',3X,'O2+N2D',
     > 3X,'O2+N4S   N2P+O2   N2A+O    N++O2    N4S+NO   O2P+NO   O++NO'
     > ,3X,'N2D+NO   hv<1910   Lyman-a')
      IF(JPR.GT.0) WRITE(I,7) Z,NNO,P1/L1,(PR(K),K=1,5)
     > ,(LR(K)*NNO,K=1,6)
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C::::::::::::::::::::::::::::::: N(4S):::::::::::::::::::::::::::::::::::::::
      SUBROUTINE CN4S(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1,L1,N4S,DISN4S
     >   ,N2D,N2P,OPLS,N2PLS,UVDISN,NOP,NPLUS,NNO,O2P,PDNOSR,VCON)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=DISN4S
      PR(2)=RTS(15)*ON*N2D
      PR(3)=RTS(8)*NE*N2D
      PR(4)=VCON*RTS(3)*OPLS*N2N
      PR(5)=RTS(53)*RTS(11)*N2PLS*NE
      PR(6)=RTS(62)*UVDISN
      PR(7)=NOP*NE*RTS(49)
      PR(8)=N2D*RTS(61)
      PR(9)=N2P*RTS(58)
      PR(10)=RTS(25)*NPLUS*O2N
      PR(11)=PDNOSR*NNO
      PR(12)=NPLUS*NNO*RTS(81)      !..Fox
      LR(1)=RTS(7)*O2N
      LR(2)=RTS(9)*NNO
      LR(3)=RTS(21)*O2P
      LR(4)=RTS(79)*N2PLS     !..Fox
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)+PR(9)+PR(10)
     >     +PR(11)+PR(12)
      L1=LR(1)+LR(2)+LR(3)+LR(4)
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,193)
 193   FORMAT(/2X,'N(4S)',38X,'PRODUCTION',46X,':',7X,'LOSS RATES'/
     > ,3X,'ALT',2X,'[N4S]',2X,'hv->N+'
     > ,3X,'O+N2D',2X,'e+N2D',3X,'O++N2',3X,'N2++e',4X,'hv->2N'
     > ,2X,'NO++e',2X,'N(2D)',4X,'N(2P)   N+&X    hv+NO    +O2  '
     > ,2X,' +NO  ',2X,'+O2+ & N2+')
      PR(10)=PR(10)+PR(12)  !... for printing fit
      LR(3)=LR(3)+LR(4)     !... for printing fit
      IF(JPR.GT.0) WRITE(I,7) Z,N4S,(PR(K),K=1,11),(LR(K)*N4S,K=1,3)
      RETURN
 7    FORMAT(F6.1,1P,22E8.1)
      END
C::::::::::::::::::::::::::::::: CN2PLS :::::::::::::::::::::::::::::::
C..... Simplified chemistry of N2+.  PUN2P* = production of N2+ by euv 
C..... in the (X,A,B states). PEN2P* same for p.e.s (X,A,B states)
      SUBROUTINE CN2PLS(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,N2PLS,PUN2PX,
     >  PUN2PA,PUN2PB,PEN2PX,PEN2PA,PEN2PB,OP2D,OP2P,HEPLUS,NPLUS,
     >  NNO,N4S)
      IMPLICIT REAL(A-H,L,N-Z)
      REAL PUN2PX,PUN2PA,PUN2PB,PEN2PX,PEN2PA,PEN2PB
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=PUN2PX
      PR(2)=PUN2PA
      PR(3)=PUN2PB
      PR(4)=PEN2PX
      PR(5)=PEN2PA
      PR(6)=PEN2PB
      PR(7)=RTS(19)*OP2D*N2N
      PR(8)=RTS(20)*OP2P*N2N
      PR(9)=RTS(44)*HEPLUS*N2N
      PR(10)=RTS(82)*NPLUS*NNO   !..Fox
      LR(1)=RTS(10)*ON
      LR(2)=RTS(11)*NE
      LR(3)=RTS(17)*O2N
      LR(4)=RTS(99)*ON
      LR(5)=RTS(79)*N4S     !..Fox
      LR(6)=RTS(80)*NNO     !..Fox
      N2PLS=
     >  (PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)+PR(9)+PR(10))/
     >  (LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,95)
 95   FORMAT(/2X,'N2+',29X,'PRODUCTION',45X,':',12X,'LOSS RATES'/
     > ,3X,'ALT  [N2+]  EUV-X   EUV-A    EUV-B   PE-X'
     > ,5X,'PE-A    PE-B  O+2D+N2  O+2P+N2  He++N2  O+N2+'
     > ,2X,'e+N2+  O2+N2+  N2++O  Other')
      PR(9)=PR(9)+PR(10)     !.. for printing fit
      LR(5)=LR(5)+LR(6)      !.. for printing fit
      IF(JPR.GT.0) WRITE(I,7) Z,N2PLS,(PR(K),K=1,9)
     > ,(LR(K)*N2PLS,K=1,5)
      RETURN
 7    FORMAT(F6.1,1P,22E8.1)
      END
C:::::::::::::::::::::::::::::: NO+ ::::::::::::::::::::::::::::::::::
      SUBROUTINE CNOP(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1,NOP,OPLS
     >  ,N2PLS,O2P,N4S,NNO,NPLUS,N2P,PLYNOP,VCON,N2D,OP2D)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=VCON*RTS(3)*N2N*OPLS
      PR(2)=N2PLS*ON*RTS(10)
      PR(3)=O2P*N4S*RTS(21)
      PR(4)=O2P*NNO*RTS(23)
      !.. N+ + O2 -> O2+ + N(2D,4S) or NO+ + O(1S)
      PR(5)=(RTS(30)+RTS(66)+RTS(59))*NPLUS*O2N
      PR(6)=RTS(37)*N2P*ON
      PR(7)=RTS(24)*OPLS*NNO
      PR(8)=PLYNOP*NNO
      PR(9)=O2P*N2D*RTS(77)         !..Fox
      PR(10)=N2PLS*NNO*RTS(80)      !..Fox
      PR(11)=NPLUS*NNO*RTS(81)      !..Fox
      PR(12)=RTS(83)*NNO*OP2D       !..Fox
      PR(13)=OP2D*RTS(90)*N2N        !.. -> NO+ + N, Li et al. [1997] 
      LR(1)=NE*RTS(5)
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)
     >    +PR(9)+PR(10)+PR(11)+PR(12)+PR(13)
      NOP=P1/LR(1)

      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,96)
 96   FORMAT(/2X,'NO+',31X,'PRODUCTION',48X,':',2X,'LOSS RATES'/
     > ,3X,'ALT',3X,'[NO+]',4X,'O++N2',3X,'N2++O',3X,'O2++N4S'
     > ,3X,'O2++NO',3X,'N++O2',4X,'N2P+O',3X,'O++NO   hv+NO'
     > ,5X,'O2++N2D   N2++NO   N++NO   OP2D+NO   OP2D+N2  NO++e')
      !PR(9)=PR(9)+PR(10)+PR(11)+PR(12)+PR(13)
      IF(JPR.GT.0) WRITE(I,7) Z,NOP,(PR(K),K=1,13),LR(1)*NOP
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C::::::::::::::::::::::::::::::: O2+ :::::::::::::::::::::::::::::::::::::::
      SUBROUTINE CO2P(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1
     > ,O2P,TPROD5,OPLS,OP2D,N2PLS,NPLUS,N4S,NNO,OP2P)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=TPROD5     !.. TPROD5=euv @ p.e. production
      PR(2)=RTS(4)*O2N*OPLS
      PR(3)=RTS(43)*OP2D*O2N
      PR(4)=RTS(17)*O2N*N2PLS
      PR(5)=RTS(25)*NPLUS*O2N
      PR(6)=RTS(86)*OP2P*O2N           !.. Fox
      PR(7)=RTS(65)*NPLUS*O2N
      LR(1)=RTS(6)*NE
      LR(2)=RTS(21)*N4S
      LR(3)=RTS(23)*NNO
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)
      O2P=P1/(LR(1)+LR(2)+LR(3))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,97)
  97  FORMAT(/2X,'O2+',22X,'PRODUCTION',24X,':',12X,'LOSS RATES'
     > /,3X,'ALT',3X,'[O2+]',3X,'hv+O2',3X,'O++O2',3X,'O+(2D)+O2'
     >  ,4X,'N2++O2   N++O2   O+(2P)+O2  O2++e   O2++N   O2++NO')
      IF(JPR.GT.0) WRITE(I,7) Z,O2P,(PR(K),K=1,7),(LR(K)*O2P,K=1,3)
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C::::::::::::::::::::::::::::::::: O+(4S) :::::::::::::::::::::::::::::::::::::
      SUBROUTINE COP4S(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,OPLS,TPROD1,OP2D
     >  ,OP2P,PEPION,PDISOP,N2PLS,N2D,NNO,VCON,HEPLUS)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=TPROD1     !.. PR(1)= euv production of o+(4s)
      PR(2)=OP2D*NE*RTS(12)
      PR(3)=OP2P*ON*RTS(26)
      PR(4)=PEPION
      PR(5)=PDISOP
      PR(6)=RTS(99)*N2PLS*ON
      PR(7)=OP2P*NE*RTS(14)
      PR(8)=OP2P*0.047
      PR(9)=RTS(28)*ON*OP2D
      PR(10)=RTS(85)*OP2P*O2N              !.. Fox
      PR(11)=HEPLUS*O2N*(RTS(91)+RTS(93))  !..Fox
      PR(12)=RTS(95)*NNO*HEPLUS            !..Fox
      PR(13)=0.0
      LR(1)=N2N*VCON*RTS(3)
      LR(2)=O2N*RTS(4)
      LR(3)=NNO*RTS(24)
      LR(4)=N2D*RTS(29)    !.... small loss?? ..Fox
      !..LR(4)=(LR(1)+LR(2)+LR(3))  !.. total loss for printing
      PR(10)=PR(10)+PR(11)+PR(12)+PR(13)
      PRTOT=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)+PR(9)+PR(10)
      LRTOT=LR(1)+LR(2)+LR(3)+LR(4)
      OPLS=PRTOT/LRTOT
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,98)
 98   FORMAT(/2X,'O+',41X,'PRODUCTION',39X,':',10X,'LOSS RATES'/
     >  ,' ALT    [O+]   hv+O  O+(2D)+e O+(2P)+O   e*+O  O2-diss  '
     >  ,'N2++O  O+(2P)+e O+(2P) O+O+(2D)   Other  +N2     +O2    '
     >  ,'+NO   +N2D')
      IF(JPR.GT.0) WRITE(I,7) Z,OPLS,(PR(K),K=1,10)
     > ,(LR(K)*OPLS,K=1,4)
      RETURN
 7    FORMAT(F6.1,1P,22E8.1)
      END
C::::::::::::::::::::::::::::::::::: O+(2D) :::::::::::::::::::::::::::::::::::
      SUBROUTINE COP2D(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,OP2D,TPROD2,OP2P
     > ,HEPLUS,N4S,NNO,PSEC)
      IMPLICIT REAL (A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=TPROD2                ! EUV  prod
      PR(2)=OP2P*NE*RTS(13)
      PR(3)=OP2P*0.171
      PR(4)=HEPLUS*O2N*RTS(76)     !..Fox
      PR(5)=PSEC
      LR(1)=RTS(19)*N2N
      LR(2)=7.7E-5                 !.. radiation at 3726 and 3729 A
      LR(3)=NE*RTS(12)
      LR(4)=ON*RTS(28)
      LR(5)=RTS(43)*O2N
      LR(6)=RTS(83)*NNO      !..Fox
      LR(7)=RTS(84)*N4S      !..Fox
      LR(8)=RTS(90)*N2N      !.. -> NO+ + N, Li et al. [1997] 
      OP2D=(PR(1)+PR(2)+PR(3)+PR(4)+PR(5))/
     >   (LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,99)
 99   FORMAT(/2X,'O+(2D)',13X,'PRODUCTION',27X,':',18X,'LOSS RATES'/
     > ,3X,'ALT',3X,'[O+2D]',3X,'hv+O',4X,'e*+O',4X,'O+2P+e',3X,
     > 'O+2P>hv',2X,'He++O2     +N2    E3726_29    +e       +O',
     > '      +O2      +NO     +N  +N2>NO+')
      IF(JPR.GT.0) WRITE(I,7) Z,OP2D,PR(1),PR(5),PR(2),PR(3),PR(4)
     > ,(LR(K)*OP2D,K=1,8)
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C::::::::::::::::::::::::::::::::::: O+(2P) :::::::::::::::::::::::::::::::::::
      SUBROUTINE COP2P(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,OP2P,TPROD3,PSEC
     > ,HEPLUS,N4S,NNO,TE)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=0.0
      IF(TPROD3.GE.PSEC) PR(1)=TPROD3-PSEC
      PR(2)=PSEC
      PR(3)=HEPLUS*O2N*RTS(92)        !..Fox
      LR(1)=RTS(26)*ON
      LR(2)=RTS(20)*N2N
      LR(3)=RTS(13)*NE
      LR(4)=0.218
      LR(5)=RTS(14)*NE
      LR(6)=(RTS(85)+RTS(86))*O2N    !..Fox
      LR(7)=RTS(87)*N4S              !..Fox
      LR(8)=RTS(88)*NNO              !..Fox
      OP2P=(TPROD3+PR(3))
     >  /(LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)+LR(8))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,100)
 100   FORMAT(/2X,' O+(2P)',6X,'PRODUCTION',10X,':',12X,'LOSS RATES'/
     >,3X,'ALT   [O+2P]    hv+O     e*+O  He++O2      +O',7X,'+N2'
     > ,6x,'+e       RAD      +e      +O2      +N4S     +NO'
     > ,6x,'OX       N2        e      Te       E7320')
      IF(JPR.GT.0) WRITE(I,7) Z,OP2P,PR(1),PR(2),PR(3),
     >  (LR(K)*OP2P,K=1,8),ON,N2N,NE,TE,OP2P*0.218*0.781
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C:::::::::::::::::::::::::::::::::: N+ ::::::::::::::::::::::::::::::::::::
      SUBROUTINE CNPLS(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,DISNP
     > ,NPLUS,OPLS,N2D,OP2P,HEPLUS,PHOTN,O2P,N4S,OP2D,N2PLS,NNO)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=DISNP
      PR(2)=RTS(29)*OPLS*N2D
      PR(3)=0
      PR(4)=RTS(45)*HEPLUS*N2N
      PR(5)=PHOTN
      PR(6)=O2P*N2D*RTS(78)          !..Fox
      PR(7)=N2PLS*N4S*RTS(79)        !..Fox
      PR(8)=OP2D*N4S*RTS(84)         !..Fox
      PR(9)=RTS(94)*NNO*HEPLUS       !..Fox
      LR(1)=RTS(30)*O2N              !..Fox
      LR(2)=RTS(25)*O2N              !..Fox
      LR(3)=RTS(22)*O2N              !..Fox
      LR(4)=RTS(65)*O2N              !..Fox
      LR(5)=RTS(66)*O2N              !..Fox
      LR(6)=RTS(31)*ON               !..Fox

      CNPLUS=0.0
      IF(LR(1)+LR(2)+LR(3).GT.0.0)
     >CNPLUS=(PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)+PR(9))
     >   /(LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6))
      NPLUS=CNPLUS

      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,101)
 101   FORMAT(/2X,'N+',20X,'PRODUCTION',71X,':',8X,'LOSS RATES'/
     > ,4X,'ALT   [N+]   [N+]c     hv+N2   O++N2D  O+2P+N2',3X
     > ,'He++N2',3X,' hv+N   O2++N2D  N2++N4S O+(2D)+N4S  He++NO'
     > ,3X,'N++O2    N++O2    N++O2    N++O2    N++O2    N++O')
      IF(JPR.GT.0) WRITE(I,7) Z,NPLUS,CNPLUS
     > ,(PR(K),K=1,9),(LR(K)*NPLUS,K=1,6)
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C::::::::::::::::::::::::::::::::: N2(A3sigma+LBH) :::::::::::::::::::::::::::::::::::::
      SUBROUTINE CN2A(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE
     > ,N2A,P3X1,P3X2,P3X3,P3X4)
      IMPLICIT REAL(A-H,L,N-Z)
      REAL P3X1,P3X2,P3X3,P3X4
      DIMENSION RTS(99),LR(22),PR(22)
      !... pr(1,2,3)= electron impact excitation of n2(a,b,c) states
      PR(1)=P3X1
      PR(2)=P3X2
      PR(3)=P3X3
      LR(1)=RTS(36)*ON
      LR(2)=RTS(27)*ON
      LR(3)=0.57
      N2A=(PR(1)+PR(2)+PR(3))/(LR(1)+LR(2)+LR(3))
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,102)
 102   FORMAT(/2X,'N2(A)',12X,'PRODUCTION',13X,':',5X,'LOSS RATES'
     > ,3X,':  Total LBH'
     > /,3X,'ALT',3X,'N2(A)',3X,'e*->N2A',3X,'e*->N2B',3X,'e*->N2C',2X
     >  ,'N2A>O1S',2X,'N2A>NO',2X,'RAD',5X,'LBH')
      IF(JPR.GT.0) WRITE(I,7) Z,N2A,(PR(K),K=1,3),(LR(K)*N2A,K=1,3)
     > ,P3X4
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C:::::::::::::::::::::::::::::::::: N(2P) ::::::::::::::::::::::::::::::::::::
C..... The rates are from Zipf et al JGR, 1980 p687
C..... 21-AUG-1992. Added N2+ recombination source
      SUBROUTINE CN2P(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1,L1
     > ,N2P,P3X7,UVDISN,O2P,NNO,N2PLUS)
      IMPLICIT REAL(A-H,L,N-Z)
      DIMENSION RTS(99),LR(22),PR(22)
      PR(1)=P3X7
      PR(2)=RTS(64)*UVDISN
      PR(3)=RTS(73)*RTS(11)*N2PLUS*NE
      LR(1)=RTS(37)*ON
      LR(2)=RTS(38)*O2N
      LR(3)=RTS(39)*O2P
      LR(4)=RTS(40)*NNO
      LR(5)=RTS(57)
      LR(6)=RTS(58)
      LR(7)=(RTS(96)+RTS(97))*NE    !..Fox
      P1=PR(1)+PR(2)+PR(3)
      L1=LR(1)+LR(2)+LR(3)+LR(4)+LR(5)+LR(6)+LR(7)
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,103)
 103   FORMAT(/2X,'N(2P)',9X,'PRODUCTION',17X,':',20X,'LOSS RATES'/
     > ,3X,'ALT',3X,'[N2P]',3X,'e+N2',5X,'hv+N2',3X,'e+N2+',6X,'+O  '
     > ,3X,'+O2      +O2+      +NO       +2D     +4S      +e')
      IF(JPR.GT.0) WRITE(I,7) Z,N2P,(PR(K),K=1,3),(LR(K)*N2P,K=1,7)
      RETURN
 7    FORMAT(F6.1,1P,22E9.2)
      END
C:::::::::::::::::::::::::::::: NO+(v) ::::::::::::::::::::::::::::::::::
C...... This routine calculates the vibrational distribution of no+
C...... Uses AFRL report by Winick et al. AFGL-TR-87-0334, Environmental
C...... Research papers, NO. 991, "An infrared spectral radiance code 
C...... for the auroral thermosphere (AARC)
C...... Written by P. Richards in February 2004
      SUBROUTINE CNOPV(JPR,I,JPT,Z,RTS,ON,O2N,N2N,NE,P1,NOP,OPLS
     >  ,N2PLS,O2P,N4S,NNO,NPLUS,N2P,PLYNOP,VCON,N2D,OP2D)
      IMPLICIT NONE
      INTEGER JPR,I,JPT,INV,IJ,IV,K
      PARAMETER (INV=20)            !.. NO+(v) array dimensions
      REAL Z,RTS(99),   !.. Altitude, rate coefficients
     >  ON,O2N,N2N,NE,              !.. O, N2, O2, electron densities
     >  P1,                         !.. total source output for finding [e] 
     >  NOP,OPLS,N2PLS,O2P,         !.. NO+, )+,N2+,O2+ densities
     >  N4S,NNO,NPLUS,N2P,          !.. N(4S), NO, N+, N(2P) densities
     >  PLYNOP,VCON,                !.. Lyman-a source, N2(v) rate factor
     >  N2D,OP2D,                   !.. N(2D), O+(2D) densities                
     >  NOPV(INV),NOPTOT,           !.. NO+(v) densities and total NO+ density
     >  LR(22),PR(22),              !.. storage for NO+ sources and sinks
     >  EINSCO1(INV),EINSCO2(INV),  !.. Einstein coeffs for delv=1,2
     >  LRV(INV),                   !.. NO+(v) + e rate factors
     >  PNOPV,LNOPV,                !.. Sources and sinks of NO+(v)
     >  PCASC,LRAD,                 !.. Temp total cascade source, sink
     >  K_N2_Q,P_N2_Q,L_N2_Q        !.. N2 queching rate :- coeff, source, sink

      !.. Fractions of each source going to each vib. level. Assume
      !.. N+ + O2 fractions for each source. NEED UPDATE
      REAL PRV1(INV),PRV2(INV),PRV3(INV),PRV4(INV),PRV5(INV),PRV6(INV),
     >  PRV7(INV),PRV8(INV),PRV9(INV),PRV10(INV),PRV11(INV),PRV12(INV)
      DATA PRV1/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV2/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV3/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV4/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV5/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV6/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV7/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV8/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV9/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV10/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV11/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/
      DATA PRV12/0,.05,.07,.09,.11,.13,.14,.17,.07,.01,.02,.06,.08,7*0/

      DATA K_N2_Q/7.0E-12/  !.. Quenching rate coeff. by N2
      DATA EINSCO1/0.0,10.9,20.2,28.4,35.5,41.5,46.6,50.9,54.3,57.0,
     >     58.9,60.2,60.8,60.7,59.9,5*60.0/  !.. Einstein coeff delv=1
      DATA EINSCO2/0.0,0.0,.697,1.93,3.61,5.74,8.24,11.1,14.2,17.7,
     >     21.3,25.1,29.0,33.2,37.4,5*40.0/ !.. Einstein coeff delv=1
      !.. rate factors for NO+(v)+e -> N + O. Sheehan and St-Maurice 2004
      DATA LRV/1.0,19*0.3333/    
      
      !... Evaluate total production and loss rates
      PR(1)=VCON*RTS(3)*N2N*OPLS        !.. N2 + O+
      PR(2)=N2PLS*ON*RTS(10)            !.. N2+ + O
      PR(3)=O2P*N4S*RTS(21)             !.. O2+ + N(4S)
      PR(4)=O2P*NNO*RTS(23)             !.. O2+ + NO
      !.. N+ + O2 -> O2+ + N(2D,4S) or NO+ + O(1S)
      PR(5)=(RTS(30)+RTS(66)+RTS(59))*NPLUS*O2N
      PR(6)=RTS(37)*N2P*ON              !.. N2+ + O
      PR(7)=RTS(24)*OPLS*NNO            !.. O+ + NO
      PR(8)=PLYNOP*NNO                  !.. Lyman-a + NO
      PR(9)=O2P*N2D*RTS(77)             !.. Fox: O2+ + N(2D)
      PR(10)=N2PLS*NNO*RTS(80)          !.. Fox: N2+ + NO
      PR(11)=NPLUS*NNO*RTS(81)          !.. Fox: N+ + NO
      PR(12)=RTS(83)*NNO*OP2D           !.. Fox: O+(2D) + NO
      PR(13)=OP2D*RTS(90)*N2N           !.. -> NO+ + N, Li et al. [1997] 
      LR(1)=NE*RTS(5)                   !.. NO+ + e

      !..Total source term used in main program to calculate [e]
      P1=PR(1)+PR(2)+PR(3)+PR(4)+PR(5)+PR(6)+PR(7)+PR(8)
     >    +PR(9)+PR(10)+PR(11)+PR(12)+PR(13)
      NOP=P1/LR(1)         !.. NO+ density

      DO IJ=1,INV
        NOPV(IJ)=0.0
      ENDDO
      NOPTOT=0.0

      !.. loop down evaluating the vibrational populations. Must start 
      !.. less than INV-1 because need cascade from v+2
      DO IV=INV-4,1,-1
        !... chemical production for v=IV = total source * fraction 
        PNOPV=PR(1)*PRV1(IV)+PR(2)*PRV2(IV)+PR(3)*PRV3(IV)+
     >    PR(4)*PRV4(IV)+PR(5)*PRV5(IV)+PR(6)*PRV6(IV)+
     >    PR(7)*PRV7(IV)+PR(8)*PRV8(IV)+PR(9)*PRV9(IV)+
     >    PR(10)*PRV10(IV)+PR(11)*PRV11(IV)+PR(12)*PRV12(IV)

        !.. cascade production from v+1, v+2
        PCASC=NOPV(IV+1)*EINSCO1(IV+1)+NOPV(IV+2)*EINSCO2(IV+2)
        LRAD=EINSCO1(IV)+EINSCO2(IV)   !.. total radiative loss

        L_N2_Q=K_N2_Q*N2N              !.. sink of quanta by N2 quenching
        IF(IV.EQ.1) L_N2_Q=0.0

        P_N2_Q=K_N2_Q*N2N*NOPV(IV+1)   !.. source of quanta by N2 quenching
        LNOPV=LR(1)*LRV(IV)            !.. recombination rate for level IV

        !.. evaluate NO+(iV) population
        NOPV(IV)=(PNOPV+PCASC+P_N2_Q)/(LNOPV+LRAD+L_N2_Q)
        NOPTOT=NOPTOT+NOPV(IV)  !.. total NO+ concentration

          !... diagnostic print. Set alt range to invoke
          IF(JPR.GT.0.AND.Z.GE.0.AND.Z.LT.10)
     >     WRITE(6,'(F10.1,I7,1P,22E10.2)') 
     >     Z,IV,PNOPV,P1,PCASC,P_N2_Q,LNOPV,LRAD,L_N2_Q,
     >     NOPV(IV),NOP,NOPTOT
      ENDDO

      PR(9)=PR(9)+PR(10)+PR(11)+PR(12)         !.. for printing only
      IF(JPT.EQ.1.AND.JPR.GT.0) WRITE(I,96)
 96   FORMAT(/5X,'NO+',34X,'PRODUCTION',69X,':LOSS RATE'/
     > ,3X,'ALT NO+(v=0) NO+(v=1) NO+(v=2)  NO+(v=3) O++N2    N2++O',
     >  4X,'O2++N4S   O2++NO   N++O2    N2P+O   O++NO   hv+NO',
     >  4X,'Other_P   NO++e')
      IF(JPR.GT.0) WRITE(I,'(F6.1,1P,22E9.2)')  Z,
     >   (NOPV(K),K=1,4),(PR(K),K=1,9),LR(1)*NOP

      RETURN
      END
C
C
C.................................... RATES.FOR ................. 
C.... This is the reaction rate subroutine for the FLIP model. It takes
C.... temperatures as input and puts the rates in array RTS. It includes
C.... reaction rates, efficiencies for various products, and Einstein
C.... coefficients. For a complete set of references see Fox and Sung JGR
C.... 2001, page 21,305. Rates different from Fox and Sung indicated by PGR
      SUBROUTINE RATS(J,TE,TI,TN,RTS)
      IMPLICIT REAL(A-H,L,N-Z)
      REAL TE,TI,TN,RTS,T13,TOT_NP_O2_RATE
      DIMENSION RTS(99)

      !.. zero out array
      DO 898 ITJS=1,99
 898  RTS(ITJS)=0.0

      !.. O + H+ -> O+ + H      Fox and Sung [2001]
      ZED=1+0.6*EXP(-228.0/TN)+0.2*EXP(-326.0/TN)
      RTS(1)=(8*6.4E-10/(9*ZED))
     >     *(EXP(-232.1/TI)+0.6*EXP(-228.0/TI)+0.2*EXP(-326.0/TI))

      !.. O+ + H -> O + H+    Anicich et al. [1993]
      RTS(2)=6.4E-10

      !.. O+ + N2 --> NO+ + N,   Hierl et al.[1997] 
      !.. The Hierl et al. [1997] lab rate is contaminated by N2(v) 
      !.. for T > 1300K. Therefore, the Hierl et al. rate is not really 
      !.. appropriate  in the ionosphere. The IDC model uses the Hierl et  
      !.. al. rate because it does not solve for N2(v). The FLIP model 
      !.. solves for N2(v) and uses the St. Maurice and Torr rate (JGR,1978,p969)
      IF(TI.LE.1000) RTS(3)=1.2E-12*(300/TI)**0.45     !.. Hierl et al.[1997] 
      IF(TI.GT.1000) RTS(3)=7.0E-13*(TI/1000)**2.12    !.. Hierl et al.[1997] 

      !.. O+ + O2 -> O2+ + O,   Lindinger et al. [1974] 
      !.. Hierl et al. lists different rates. Hierl et al. [1997] not 
      !.. used above 1600 because rates are contaminated by O2(v) for 
      !.. T > 1000K. We don't know the vibrational state in the 
      !.. thermosphere. This fit was done by PGR May 2009. It is similar 
      !.. to Fox and Sung but does not increase sharply above 1000K.
      IF(TI.LE.1600) RTS(4)=1.6E-11*(300/TI)**0.52
      IF(TI.GT.1600) RTS(4)=6.7E-12*(TI/1600)**0.6 

      !.. NO+ + e -> N + O    Walls and Dunn [1974)
      !.. Vejby-Christensen et al [1998] gives 4.0E-7*(300/TE)**0.5
      !.. Torr and Torr [1979] gives 4.3E-7*(300/TE)**0.83(-0.16,+.08)
      !.. Sheehan and St. Maurice gives 3.5E-7*(300/TE)**0.65
      RTS(5)=4.0E-7*(300/TE)**0.85

      !.. O2+ + e -> O + O   Mehr and Biondi (1969)
      IF(TE.LE.1200) RTS(6)=1.953E-7*(300/TE)**0.70
      IF(TE.GT.1200) RTS(6)=7.389E-8*(1200/TE)**0.56

      !..  O2 + N(4S)-> NO + O           Baulch et al.[1994]
      RTS(7)=1.5E-14*TN*EXP(-3270.0/TN)

      !..   N(2D) + e -> N(4S) + e     Berrington and Burke [1981]
      RTS(8)=3.86E-10*(TE/300.)**0.81

      !..   NO + N(4S) -> N2 + O      Lee et al. [1978]
      RTS(9)=3.4E-11 

      !.. N2+ + O -> NO+ + N   Scott et al.[1999]
      IF(TI.LE.1500) RTS(10)= 1.33E-10*(300/TI)**0.44
      IF(TI.GT.1500) RTS(10)= 6.55E-11*(1500/TI)**(-0.2)

      !.. N2+ + e -> N + N  Mehr and Biondi (1969)
      RTS(11)=2.2E-7*(300/TE)**0.39    !.. Zipf (1980)

      !.. O+(2D) + e -> O+(4S) + e   McLaughlin and Bell (1998)
      !.. Henry [1969] gives 7.8E-8*(300/TE)**0.5
      RTS(12)=6.03E-8*(300/TE)**0.5

      !.. O+(2P) + e ->  O+(2D) + e   McLaughlin and Bell (1998)
      !.. RTS(13)+RTS(14) agrees with Walker et al (1975) and 
      !.. Chang et al (1993)
      RTS(13)=1.84E-7*(300/TE)**0.5

      !.. O+(2P) + e -> O+(4S) + e  McLaughlin and Bell (1998)
      RTS(14)=3.03E-8*(300/TE)**0.5

      !.. N(2D) + O ->  N(4S) + O  Fell et al.[1990]. Lin and Kaufman[1971]
      RTS(15)=6.9E-13

      !.. N(2D) + O2 -> NO + O  Herron[1999]. Shihira et al.[1994]
      RTS(16)=9.7E-12*EXP(-185.0/TN)

      !.. N2+ + O2 -> O2+ + N2   Scott et al.[1999]
      IF(TI.LT.1000) THEN
        RTS(17)=5.1E-11*(300/TI)**1.16
      ELSEIF(TI.LE.2000) THEN
        RTS(17)=1.26E-11*(TI/1000)**0.67
      ELSE
        RTS(17)=2.39E-11
      ENDIF

      !.. thermal electron excitation of O(1D); Rees et al 1967 pss, p1097 .....
      RTS(18)=1.1E-10*SQRT(TE)*EXP(-2.27E+4/TE)*(0.406
     >  +0.357E-4*TE-(0.333+0.183E-4*TE)*EXP(-1.37E4/TE)
     >  -(0.456+0.174E-4*TE)*EXP(-2.97E4/TE))

      !.. N2 + O+(2D) -> N2+ + O   
      !..RTS(19)=8.0E-10                   !.. Johnson and Biondi
      RTS(19)=1.50E-10*(300/Ti)**(-0.55)   !.. Li et al by PGR
      
      !.. N2 + O+(2P) -> N2+ + 0    Fox 
      !.. RTS(20)=6.2E-10*EXP(-340/TI)   !.. Li et al from Fox wrong
      RTS(20)=2.0E-10*(300/Ti)**(-0.55)    !.. Li et al by PGR

      !.. O2+ + N(4S) -> NO+ + 0   Scott et al.[1999]
      RTS(21)=1.0E-10

      !.. N+ + O2 -> O+ + NO 
      !.. Torr and Torr gives 6.0E-10 for total N+ + O2 reaction rate
      !.. Dotan et al [1997] from Fox and Sung gives
      !IF(TI.LE.1000) TOT_NP_O2_RATE=2.02E-10*(300/TI)**(-0.45)
      !IF(TI.GT.1000) TOT_NP_O2_RATE=3.49E-10
      !.. does not seem to be correct. Probably vibrationally excited O2
      !.. Branching ratios for N+ + O2 from O'Keefe et al J. Chem. Phys. 1986
      !.. NO+ +O(3P) = .09, NO+ + O(1D) = .36, O2+ + N(4S) = 0.35, 
      !.. O2+ + N(2D) = 0.15, O+(4S) + NO = .05
      TOT_NP_O2_RATE=6.0E-10                !.. Total N+ + O2 rate
      RTS(22)=0.05*TOT_NP_O2_RATE      

      !.. O2+ + NO -> NO+ + O2 Midey and Viggiano [1999]
      RTS(23)=4.5E-10 * 1.0000

      !.. O+ + NO -> O + NO+   Dotan and Viggiano [1999]
      IF(TI.LE.300) RTS(24)=7.0E-13*(300/TI)**0.66
      IF(TI.GT.300) RTS(24)=7.0E-13*(TI/300)**0.87

      !.. N+ + O2 -> O2+ + N(4S) 
      RTS(25)=0.35*TOT_NP_O2_RATE 

      !.. O+(2P) + O -> O+(4S) + O 
      !..RTS(26)=5.2E-10  !.. Fox appears to be wrong  
      !.. (Chang et al., JGR 1993) c.f. 5.2E-11  (Rusch)    
      RTS(26)=4.0E-10

      !.. N2(A3sig) + O -> NO + N(2D) 
      RTS(27)=2.0E-11       !..see Campbell et al. 2006
      RTS(27)=0.000000      !.. Torr and Torr value

      !.. O+(2D) + O ->  O+(4S) + O  Torr and Torr [1980]
      RTS(28)=1.0E-11

      !.. O+ + N(2D) -> O + N+  Constantinides et al.[1979].Bates[1989]
      RTS(29)=1.3E-10

      !.. O2 + N+ -> O(3P) + NO+
      !.. Branching ratio from O'Keefe et al J. Chem. Phys. 1968
      RTS(30)=0.09*TOT_NP_O2_RATE 

      !.. O + N+ -> O+ + N   Constantinides et al.[1979].Bates[1989]
      RTS(31)=2.2E-12

      !.. Efficiency for   N2+ + e -> N(2D) + N(2D)
      RTS(32)=1.46

      !.. N2 + O(1D) -> O + NO
      RTS(33)=1.8E-11*EXP(107.0/TN)

      !.. O2 + O(1D) -> O + O2
      RTS(34)=3.2E-11*EXP(67/TN)

      !.. O2 + N(4S) -> O(1S) + NO. Kopp et al. 1977, JGR, p4715
      RTS(35)=2.5E-11


      !.. N2(A3sig) + O -> O(1S) + N2
      RTS(36)=2.5E-11*EXP(TN/298)**0.55    !..  see Campbell et al. 2006
      RTS(36)=2.0E-11                      ! .. Torr et al.

      !.. N(2P) + O -> products (N(2D,4S) and NO+) and O(3P,1D) 
      !.. from Piper et al 1993, J. Chem. Phys. vol 98 page 8560.
      RTS(37)=1.7E-11

      !.. N(2P) + O2 -> NO + O 
      RTS(38)=3.9E-12*EXP(-60/TN)

      !.. N(2P) quenching rates(O2+,NO) from Zipf et al jgr 1980 p687
      RTS(39)=2.2E-11
      RTS(40)=1.8E-10

      !.. N(2D) + NO -> N2 + O 
      RTS(41)=6.7E-11

      !.. efficiency N2+ + O -> N2 + O+(4S)   
      IF(TI.LE.1500) RTS(42)= 7.0E-12*(300/TI)**0.21
      IF(TI.GT.1500) RTS(42)= 4.83E-12*(1500/TI)**(-0.41)
      RTS(42)=RTS(42)/RTS(10)    !.. converts to efficiency

      !.. O+(2D) + O2 -> O2+ + O   Fox
      RTS(43)=7.0E-10

      !.. He+ + N2 -> He + N2+
      RTS(44)=5.2E-10

      !.. He+ + N2 -> He + N+
      RTS(45)=7.8E-10

      !.. O(1S)+ e -> O(1D) + e  
      RTS(46)=8.5E-9

      !.. O(1S)+ e -> O(3P) + e  
      RTS(47)=1.56E-10*(TE/300)**0.94

      !.. O(1S) + O2 -> O2 + O 
      RTS(48)=4.4E-12*EXP(-815.0/TN)

      !.. NO+ + e -> N(4S) + O
      RTS(49)=0.15 * RTS(5)

      !.. NO+ + e -> N(2D) + O
      RTS(50)=0.85 * RTS(5)

      !.. O2+ + e -> O(1D) + O
      RTS(51)=1.11 * RTS(6)

      !.. O2+ + e -> O(1S) + O
      RTS(52)=0.05 * RTS(6)

      !.. Efficiency for   N2+ + e -> N(4S) + N(2D)
      RTS(53)=0.46

      !.. O(1D) -> O + 6300 + 6364 
      RTS(54)=0.00934

      !.. O(1S) -> O(1D) + 5577
      RTS(55)= 1.06

      !.. O(1S) -> O(3P) + hv (2972) RTS(56)= 4.5E-2 !.. old value
      RTS(56)= 0.10 * RTS(55)  !.. From Slanger, Spring AGU 2005

      !.. N(2P) -> N(2D) + hv
      RTS(57)=7.9E-2

      !.. N(2P) -> N(4S) + hv
      RTS(58)=5.0E-3

      !.. N+ + O2 -> NO+ + O(1S) Langford et al., PSS, 33,1225,1985
      RTS(59)=1.0E-3*TOT_NP_O2_RATE 

      !.. Efficiency for   N2(A3sig) + O -> O(1S) + N2
      RTS(60)=0.37

      !.. N(2D) -> N(4S) + hv
      RTS(61)=1.07E-5

      !.. hv(>600A) + N2 -> N(4S) + N   branching ratio
      RTS(62)=0.5

      !.. hv(>600A) + N2 -> N(2D) + N   branching ratio
      RTS(63)=0.4

      !.. hv(>600A) + N2 -> N(2P) + N   branching ratio
      RTS(64)=0.1

      !.. N+ + O2 -> O2+ + N(2D)
      !.. Branching ratio from O'Keefe et al J. Chem. Phys. 1968
      RTS(65)=0.15*TOT_NP_O2_RATE 

      !.. N+ + O2 -> NO+ + O(1D)
      !.. Branching ratio from O'Keefe et al J. Chem. Phys. 1968
      RTS(66)=0.36*TOT_NP_O2_RATE 

      !.. hv(Scum-Runge) + O2 -> O(1S) + O   branching ratio
      RTS(67)=0.001

      !.. Effic of O2(A3,DEL) + O -> O(1S)
      RTS(68)=0.1

      !.. O(1D) + O -> O + O   Abreu et al. PSS, p1143, 1986
      RTS(69)=6.47E-12*(TN/300)**0.14

      !.. hv + N2 -> N+(5S) -> 2143 A emission yield from the 2s sigma g state  
      !.. of N2. This was taken as 0.6 the value of Cleary and Barth JGR 1987, 
      !.. p13,635 because they did not double EUV below 250 A.
      RTS(70)=0.06

      !.. hv + N2 -> N+(1D) -> 6584 A emission (guess)
      RTS(71)=0.3

      !.. hv + N2 -> N+(1S) -> 5755 A emission (guess)
      RTS(72)=0.03

      !.. efficiency of production of N(2P) from e + N2+ reaction
      RTS(73)=0.08

      !.. Efficiency for production of O(1D) from N(2D) + O2 reaction
      !.. See Link and Swaminathan, PSS, 1992, page 699
      RTS(74)=0.1    !??? check

      !.. He+ + O2 -> He + O2+
      RTS(75) = 9.2E-12

      !.. He+ + O2 -> He + O+(2D) + O(3P) 
      RTS(76) = 2.37E-10

      !.. O2+ + N(2D) -> NO+ + O
      RTS(77) = 1.8E-10

      !.. O2+ + N(2D) -> N+ + O2
      RTS(78) = 8.65E-10

      !.. N2+ + N(4S) -> N+ + N2
      RTS(79) = 1.0E-11

      !.. N2+ + NO -> NO+ + N2
      RTS(80) = 3.6E-10

      !.. N+ + NO -> N(4S) + NO+
      RTS(81) = 4.72E-10*(300/TI)**0.24

      !.. N+ + NO -> N2+ + O
      RTS(82) = 8.33E-11*(300/TI)**0.24

      !.. O+(2D) + NO -> NO+ + O
      RTS(83)=1.2E-9

      !.. O+(2D) + N -> N+ + O
      RTS(84)=1.5E-10

      !.. O+(2P) + O2 -> O+ + O2  Fox
      RTS(85)=1.3E-10

      !.. O+(2P) + O2 -> O2+ + O
      RTS(86)=1.3E-10

      !.. O+(2P) + N -> O+ + N(2D)
      RTS(87)=1.0E-11

      !.. O+(2P) + NO -> NO+ + O
      RTS(88)=1.2E-9

      !.. H+ + O2 -> O2+ + H
      RTS(89)=3.8E-9
 
      !.. O+(2D) + N2 -> NO+ + N  !.. Li et al. (1997).
      !.. From the ratio of the cross sections.
      !.. The branching ratio to O+(4S) + N2 not given by Li et al.
      RTS(90)=2.5E-11

      !.. He+ + O2 -> He + O+(4S) + O 
      RTS(91) = 2.39E-11

      !.. He+ + O2 -> He + O+(2P) + O 
      RTS(92) = 6.04E-10

      !.. He+ + O2 -> He + O+(4S) + O(1D)
      RTS(93) = 4.6E-11

      !.. He+ + NO -> He + N+ + O
      RTS(94) = 1.35E-9

      !.. He+ + NO -> He + O+ + N
      RTS(95) = 1.0E-10

      !.. N(2P) + e -> N(4S) + e
      RTS(96)=2.04E-10*(TE/300)**0.85
 
      !.. N(2P) + e -> N(2D) + e
      RTS(97)=9.5E-9

      !.. O(1D) + e -> O(3P) + e
      RTS(98)=2.87E-10*(TE/300)**0.91

      !.. N2+ + O -> O+ + N2  !.. McFarland et al.(1974)
      !.. From Fox and Sung 2001
      RTS(99)=0.07*1.0E-10*(300/Ti)**0.23

        RETURN
        END
C
C
C.........................<Photoelectron model>......................3 APR 92
C..... Calculate secondary ion production, electron heating rate and 3371 excitation rate.
      SUBROUTINE SECIPRD(ALT,SZADEG,F107,F107A,TE,TN,OXN,O2N,N2N,XNE,
     >   N2APRD)
      INTEGER I,K,IK           !-- loop control variables
      INTEGER IDIM             !.. Array dimensions
      INTEGER IMAX             !.. Maximum PE energy
      PARAMETER (IDIM=501)
      REAL PEFLUX(IDIM)        !-- PE flux
      REAL SIGIT(3)            !-- Total electron impact cross sections
      REAL SIGEX(22)           !-- Electron impact cross sections for OX
      REAL ALT                 !-- ALT = altitude (km)  { 120 -> 500 }
      REAL SZADEG              !-- solar zenith angle {0 -> 90 degrees}
      REAL F107, F107A         !-- F107 = Solar 10.7 cm flux
      REAL TE,TN               !-- electron, neutral temperatures (K)
      REAL XN(3),OXN,O2N,N2N   !-- XN, O, O2, N2, neutral densities  (cm-3)
      REAL XNE                 !-- electron density  (cm-3)
      REAL XN2D                !-- N(2D) density for N(2D) + e -> 2.5 eV
      REAL XOP2D               !-- O+(2D) density for O+(2D) + e -> 3.3 eV
      REAL SPRD(3,6)
      REAL SIGOX,SIGN2,SIGEE  !.. Total exciation cross sections for O, N2, O2
      REAL N2APRD             !.. Production of N2A
      REAL DE(IDIM),EV(IDIM)
      !.. various ionization and excitation rates by EUV and PE
      REAL EUVION,PEXCIT,PEPION,OTHPR1,OTHPR2
      COMMON/EUVPRD/EUVION(3,12),PEXCIT(3,12),PEPION(3,12),OTHPR1(6)
     >   ,OTHPR2(6)

      DATA SPRD/.4,.56,.44, .4,.28,.44, .2,.06,.10, 0.,.05,.00, 0.,.05
     >             ,.00, 0.0,0.0,0.02/
      DATA IMAX/0/              !.. Initialize IMAX Reset in FLXCAL
      
      IMAX=0   ! IMAX needs to be initialize at each call to subroutine 
      !.. Transfer neutral densities to the density array
      XN(1)=OXN
      XN(2)=O2N
      XN(3)=N2N
      N2APRD=0.0
      DO K=1,3
      DO IK=1,6
        PEPION(K,IK)=1.0E-15
        PEXCIT(K,IK)=1.0E-15
      ENDDO
      ENDDO

      !.. Cannot calculate PE if no densities
      IF(OXN.LT.1.0E5.OR.N2N.LT.1.0E5) RETURN
      IF(SZADEG.GT.105) RETURN

      !********************************************************************
      !.. Go and get the photoelectron fluxes
       XN2D=0        !.. N(2D) density for calculating N(2D) + e -> 2.5 eV
       XOP2D=0       !.. O+(2D) density for calculating O+(2D) + e -> 3.3 eV
       CALL FLXCAL(IDIM,ALT,SZADEG,
     >   TE,TN,XN,XNE,XN2D,XOP2D,PEFLUX,AFAC,IMAX,DE,EV)
      !***************************************************************

      !........ sample calculation of ion production rates. 
      DO I=1,IMAX
        E=EV(I)
        CALL TXSION(E,SIGIT)                     !.. total ion XS
        CALL SIGEXS(E,TE,XNE,SIGOX,SIGN2,SIGEE)  !.. Total excitation XS
        CALL OXSIGS(E,SIGEX,SIGOX)               !.. OX cross sections

        IF(E.LT.250) N2APRD=N2APRD+0.22*PEFLUX(I)*SIGN2*XN(3)*DE(I) !.. N2(A) prod
        PEXCIT(1,1)=PEXCIT(1,1)+PEFLUX(I)*SIGEX(1)*XN(1)*DE(I)   !.. O(1D) prod
        PEXCIT(1,2)=PEXCIT(1,2)+PEFLUX(I)*SIGEX(2)*XN(1)*DE(I)   !.. O(1S) prod

        !.. Evaluate ionization branching ratios for O+
        CALL OXRAT(E,SPRD(1,1),SPRD(1,2),SPRD(1,3))

        !.. Calculate ion production rates
        DO K=1,3
        DO IK=1,6
          PEPION(K,IK)=PEPION(K,IK)+PEFLUX(I)*SIGIT(K)*XN(K)*SPRD(K,IK)*
     >      DE(I)
        ENDDO
        ENDDO

      EP=E+17
      PEFLX=PEFLUX(I)/12.57
c      WRITE(29,'(3F8.1,1P,22E10.2)') ALT,E,12398/EP,PEFLX,
c     >     PEFLUX(I),(SIGIT(K),K=1,3),T_XS_OX(EP),2.2*T_XS_OX(EP),
c     >     T_XS_N2(EP),PEPION(1,1),PEXCIT(1,1)
      ENDDO

      RETURN
      END
C:::::::::::::::::::::::::: PHOTOELECTRON MODEL  ::::::::::::::::::::::::
C....... This subroutine evaluates the photoelectron flux using the concept
C.......  production frequencies developed by Phil Richards at Utah 
C....... State University March 1984. It supercedes the model described in
C....... JGR, p2155, 1983. Contact EAST::CSPARA::RICHARDS on SPAN network
C------- Some minor updates in April 1992 indicated by C----
C....... I would appreciate any feedback on bugs or clarity and if it 
C....... contributes substantially to a paper, I would appreciate the 
C....... appropriate acknowledgement.
C......       **************** WARNING ****************
C...... This program is constructed to produce reasonable agreement with
C...... the Atmosphere Explorer-E PES fluxes of John Doering (Lee et al.
C...... PSS 1980, page 947). It will NOT give good fluxes if the EUV 
C...... attenuation is greater than about a factor of 7 (AFAC < 0.14).
C...... The model accurately reproduces the measured fluxes very closely
C...... for the case in the test driver at 148 km SZA=53 when AFAC=0.19.
C...... You should compare the output against the Lee et al. 1980 fluxes
C...... periodically as a check. It is doubtful below 140km during the
C...... day and below 200km near sunset. Between 200km & 350km, it should
C...... be good for solar zenith angles < 90 degrees. Above 350 km there
C...... is considerable uncertainty due to neglect of transport but most
C...... models have similar uncertainties at high altitudes due to the 
C...... uncertainty in the conjugate photoelectron flux, and the pitch 
C...... angle distribution.
C
C------ ALT = altitude (km)  { 120 -> 500 }
C------ SZADEG = solar zenith angle  {0 -> 90 degrees ? }
C------ TE, TN = electron, neutral temperatures (K)
C------ XN, XNE = O, O2, N2, and electron densities  (cm-3)
C------ XN2D, XOP2D = N(2D) and O+(2D) densities for electron quenching
C------ (cm-3). You may put these to ZERO if not available.
C------ PEFLUX = photoelectron flux to be returned (eV cm2 sec)-1
C------ AFAC = the solar EUV attenuation warning flag
      SUBROUTINE FLXCAL(IDIM,ALT,SZADEG,
     >  TE,TN,XN,XNE,XN2D,XOP2D,PEFLUX,AFAC,IMAX,DE,EV)
      INTEGER RDIM,IMAX
      REAL EMAX             !.. maximum photoelectron energy.
      PARAMETER (RDIM=84)
c      REAL RJOX(RDIM),RJN2(RDIM),XN(3),COLUM(3),PEFLUX(IDIM)
c     >  ,DE(RDIM),DELTE(RDIM),EV(RDIM),EN(RDIM),UVFAC,EUV
c      COMMON/SOL/UVFAC(59),EUV
      REAL RJOX(RDIM),RJN2(RDIM),XN(3),COLUM(3),PEFLUX(IDIM)
     >  ,DE(RDIM),DELTE(RDIM),EV(RDIM),EN(RDIM),UVFAC(59),EUV
      COMMON/SOL/UVFAC,EUV

      !.. photoelectron production frequencies by 1.0E9. Renormalized below
      !-- O production frequencies
      DATA RJOX/10*19,15,18,14,10,13,9,13,9,7,11,6,26,6,31,6,5,22,4,
     >   4,5,4.04,2.56,1.9,2.28,2.12,0.96,0.24,0.14,0.14,0.1,0.1,0.1
     >   ,0.1,0.1, 10*.05,10*.04,10*.01,10*.02/

      !-- N2 production frequencies
      DATA RJN2/6*40,43,35,35,28,29,21,25,19,19,13,19,16,12,11,7,18,8,
     >  46,27,5,5,5,5,5,5.34,2.92,1.84,2.22,1.62,0.62,0.05,0.05,0.05,
     >  0.05,0.05,0.05,0.05,0.044,10*.02,10*.01,10*.01,10*.01/

      !-- PE energy grid
      DATA EN/0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,
     >  13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5,22.5,23.5,24.5,
     >  25.5,26.5,27.5,28.5,29.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,
     >  67.5,72.5,77.5,82.5,87.5,92.5,97.5,105,115,125,135,145,155,
     >  165,175,185,195,205,215,225,235,245,255,265,275,285,295,305,
     >  315,325,335,345,355,365,375,385,395,405,415,425,435,445,455,
     >  465,475,485,495/

      !-- PE energy steps
      DATA DELTE/30*1.0,14*5.0,40*10/
      DATA EMAX/286.0/          !..  Maximum PE energy

      SZA = SZADEG/57.29578   !.. convert solar zenith angle to radians

      IF(IMAX.LT.10) THEN
        !.. transfer energy grid to pass back
        DO IE=1,RDIM
          IF(EN(IE).LT.EMAX) IMAX=IE
          DE(IE)=DELTE(IE)
          EV(IE)=EN(IE)
        ENDDO
      ENDIF

      !.. 2.5eV production from electron quenching of N2D
      PN2D=XN2D*XNE*6.0E-10*SQRT(TE/300.0)
      !.. 3.3eV production from electron quenching of O+(2D)
      POP2D=XOP2D*XNE*6.6E-8*SQRT(300./TE)

      CASEL=0.0

      !.. evaluate the neutral column density  &&&&&&&&
      CALL SCOLUM(I,SZA,ALT*1.0E5,TN,XN,COLUM)

      !...... begin flux calculation loop............................
      DO 133 IE=1,IMAX
        I=IMAX+1-IE
        IF(I.LT.1) GO TO 55
        PEFLUX(I)=0.0

        !... evaluate energy of photon responsible for electron at energy EE
        EE=EV(I)
        EP=EE+17
        IF(EE.LT.22) EP=45
        IF(EE.GE.22.AND.EE.LT.28) EP=41
        IF(EE.GE.28.AND.EE.LT.38) EP=49

        !.. evaluate total photoionization cross sections for photon energy EP
        XSOXT=T_XS_OX(EP)         !.. New OX cross section
        XSO2T=2.2*T_XS_OX(EP)     !.. O2 XS is 2.2* O XS
        XSN2T=T_XS_N2(EP)         !.. New N2 cross section

        !... evaluate EUV attenuation factor AFAC
        TAU=COLUM(1)*XSOXT+COLUM(2)*XSO2T+COLUM(3)*XSN2T
        AFAC=EXP(-TAU)

        !..... low energy cascade production from O(1D) and N2* impact
        CASOX=0.0
        IF(EE.LT.10) CASOX=PEFLUX(I+2)*SIGOX*XN(1)
        CASN2=0.0
        IF(EE.LT.6) CASN2=PEFLUX(I+1)*SIGN2*XN(3)

        !..... cascade production from thermal electron degradation
        CASEL=0.0
        IF(I.LT.IMAX) CASEL=PEFLUX(I+1)*SIGEE*XNE

        !... Production from electron quenching of metastables
        EPN2D=0.0
        IF(NINT(EE).EQ.3) EPN2D=PN2D
        EPOP2D=0.0
        IF(NINT(EE).EQ.4) EPOP2D=POP2D

        !.... evaluate cross sections (must be after cascade production)
        CALL SIGEXS(EE,TE,XNE,SIGOX,SIGN2,SIGEE)

        !..... adjust EUV production rate for different period of solar cycle
        CALL FACFLX(EE,UVFAC,FFAC)

        !..... Production of pe's at energy EE, taking into account
        !..... attenuation and EUV variation, and renormalize frequencies

        PRODOX=RJOX(I)*XN(1)*AFAC*FFAC*1.0E-9 
        PRODN2=RJN2(I)*XN(3)*AFAC*FFAC*1.0E-9 

        !..... Sum all the production rates
        PROD=PRODOX+PRODN2+CASEL+CASOX+CASN2+EPN2D+EPOP2D

        !.... total loss through collisions
        RLOSS=SIGOX*XN(1)+SIGN2*XN(3)+SIGEE*XNE

        !....... evaluate photoelectron flux
        PEFLUX(I)=PROD/RLOSS
 133  CONTINUE

 55   RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE FACFLX(EE,UVFAC,FFAC)
C....... solar UVFAC factors. Correspond to the first 9 wavelengths
C....... TORR et al.[1979] GRL page 771 table 3. UVFAC(9) is for 304A
      REAL UVFAC(59)
      FFAC=(7*UVFAC(9)+UVFAC(8)+0.2*UVFAC(6))/8.2
      IF(EE.GT.30.AND.EE.LE.38) FFAC=(2*UVFAC(7)+.5*UVFAC(5))/2.5
      IF(EE.GT.38.AND.EE.LE.45) FFAC=UVFAC(4)
      IF(EE.GT.45.AND.EE.LE.66) FFAC=UVFAC(3)
      IF(EE.GT.66.AND.EE.LE.108) FFAC=UVFAC(2)
      IF(EE.GT.108) FFAC=UVFAC(1)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SIGEXS(E,TE,XNE,SIGOX,SIGN2,SIGEE)
C..... Program for evaluating the total inelastic cross sections
C
C........ loss to thermal electrons ....
      ET=8.618E-5*TE
      SIGEE=(3.37E-12/E**0.94/XNE**0.03)*((E-ET)/
     >   (E-(0.53*ET)))**2.36
C
C...... cross section for o(1d)
      SIGO1D=0.0
      IF(E.GT.1.96) SIGO1D=4E-16*(1-1.96/E)**2/E
C...... total excitation cross section for O excluding O(1D)
      IF(E.LT.25) SIGO=(0.4*E-5)*1.4E-17
      IF(E.GE.25) SIGO=7.0E-17
      IF(SIGO.LT.0.0) SIGO=0.0
C
C...... total excitation cross section for N2......
      IF(E.LT.12) SIGN2=(15.5*E-104.8)*1.7E-18
      IF(E.LT.4.0) SIGN2=5.0E-9*(1-1.4/E)**9 * (1.4/E)**16
      IF(E.GT.11.5) SIGN2=1.4E-16
      IF(SIGN2.LT.0.0) SIGN2=0.0
C
C........ total ionization cross sections from Keiffer and Dunn ....
      SIGION=0.0
      AL=ALOG10(E)
      IF(AL.LT.2.7.AND.AL.GE.1.2) SIGION=-3.6E-16*(AL-1.2)*(AL-3)
      IF(AL.GT.2.7) SIGION=1.2E-14*EXP(-AL*1.6)
      IF(E.LT.50) SIGION=1.0E-16*(0.068*E-1.06)
      IF(SIGION.LE.0.0) SIGION=0.0
C
      SIGOX=SIGO1D+SIGO+0.5*SIGION
      SIGN2=SIGN2+SIGION
      RETURN
      END
C::::::::::::::::::::::: TXSION ::::::::::::::::::::::::::::::::::
C..... total ionization cross sections for O, O2, and N2
C..... ionization cross sections keiffer and dunn ........
C..... The N2+ and O2+ cross sections were modified in April 99 to agree
C..... with the Schram et al. cross sections at high energies
      SUBROUTINE TXSION(E,SIGIT)
      DIMENSION SIGIT(3)

      !... SIGTMP is used for N2+ and O2+ at the high energies
      SIGTMP=1.0E-13*EXP(-2.303*ALOG10(E))

      !... N2+ cross section
      SIGIT(3)=0.0
      IF(E.GT.15.0) SIGIT(3)=1.42E-14*(1-9.0/E)**7.1*E**(-0.7) 
      IF(SIGTMP.LT.SIGIT(3)) SIGIT(3)=SIGTMP
      !... This correction to convert units to cm**2. Keiffer and Dunn page 10
      SIGIT(3)=0.87972*SIGIT(3)

      !... O2+ cross section
      SIGIT(2)=0.0
      IF(E.GT.12.0) SIGIT(2)=1.08E-14*(1-7.0/E)**8.6*E**(-0.65)
      IF(SIGTMP.LT.SIGIT(2)) SIGIT(2)=SIGTMP
      !... This correction to convert units to cm**2. Keiffer and Dunn page 10
      SIGIT(2)=0.87972*SIGIT(2)

      !... O+ cross section from Brook et al. J. Phys. B. Vol 11 p 3115, 1978
      SIGIT(1)=0.0
      IF(E.GT.12.0) SIGIT(1)=7.33E-15*(1-2.0/E)**34.3*E**(-0.7)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::: OXRAT ::::::::::::::::::::::::::::::::
      SUBROUTINE OXRAT(E,R4S,R2D,R2P)
C....... This subroutine returns the electron impact branching ratios
C....... for atomic oxygen from Burnett and Rountree Phys. Rev. A. 20
C....... 1979 page 1468
      R4S=1.0
      R2D=0.0
      R2P=0.0
      EV=E
      IF(E.GT.100.0) EV=100.0
      IF(EV.GT.17) R4S=-1.6E-3*EV+0.56
      IF(EV.GT.17) R2D=1.067E-3*EV+0.2933
      R2P=1-R4S-R2D
         IF(EV.LT.22) THEN
         R2P=0.0
         RTOT=R4S+R2D
         R4S=R4S/RTOT
         R2D=R2D/RTOT
         ENDIF
      RETURN
      END
C::::::::::::::::::::: T_XS_N2 :::::::::::::::::::::::::::
C.... This function calculates the N2 total photoionization
C.... cross section. P. Richards 2003-10-04
      REAL FUNCTION T_XS_N2(EP)
      IMPLICIT NONE
      REAL EP   !... photon energy
      REAL ESAVE
      DATA ESAVE/0.0/

      !.. Wavelength < 20 A, Auger ionization
      IF(EP.GE.600.0) THEN              
        T_XS_N2=0.5E-18
      !.. Wavelength < 31 A, Auger ionization
      ELSEIF(EP.GE.400.0) THEN              
        T_XS_N2=1.0E-18
      !.. Wavelength 31.62 to 23.70 A
      ELSEIF(EP.GE.392.0) THEN
        T_XS_N2=EXP(7.9864*ALOG(EP)-91.6604)
      !.. Wavelength 225 to 125 A
      ELSEIF(EP.GE.55.09) THEN
        T_XS_N2=EXP(-2.3711*ALOG(EP)-29.8142) 
      !.. Wavelength > 225 A
      ELSE
        T_XS_N2=EXP(-1.1077*ALOG(EP)-34.8787)  
      ENDIF

      !..IF(NINT(10*EP).NE.NINT(10*ESAVE)) WRITE(6,'(2F8.1,1P,2E10.2)') 
      !..> 12394.224/EP,EP, T_XS_N2/(3.39E-17*EXP(-0.0263*EP)), T_XS_N2
       ESAVE=EP

      !.. old parameterization
      !..T_XS_N2=3.39E-17*EXP(-0.0263*EP)

      RETURN
      END
C::::::::::::::::::::: T_XS_OX :::::::::::::::::::::::::::
C.... This function calculates the OX total photoionization
C.... cross section. P. Richards 2003-10-04
C.... Samson and Pareek Phys. Rev. A, 31, 1470, 1985

      REAL FUNCTION T_XS_OX(EP)
      IMPLICIT NONE
      REAL EP   !... photon energy
      REAL ESAVE
      DATA ESAVE/0.0/

      !.. NEW parameterization
      IF(EP.GE.500.0) THEN                 
        !.. Wavelength shorter than 25 A, Auger ionization
        T_XS_OX=0.5E-18
      ELSEIF(EP.GE.165.26) THEN                 
        !.. Wavelength shorter than 75 A
        T_XS_OX=EXP(-2.5209*ALOG(EP)-28.8855)
      ELSEIF(EP.GE.55.09) THEN              
        !.. Wavelength between 78 and 256.26 A
        T_XS_OX=EXP(-1.7871*ALOG(EP)-32.6335)
      ELSE
        !.. Wavelength longer than 256.26 A
        T_XS_OX=EXP(-1.3077*ALOG(EP)-34.5556)   
      ENDIF

      !..IF(NINT(10*EP).NE.NINT(10*ESAVE)) WRITE(6,'(2F8.1,1P,2E10.2)') 
      !..> 12394.224/EP,EP, T_XS_OX/(27.2E-18*EXP(-3.09E-2*EP)), T_XS_OX
       ESAVE=EP

      !.. old parameterization
      !.. T_XS_OX=27.2E-18*EXP(-3.09E-2*EP)

      RETURN
      END
C
C
C:::::::::::::::::::::: OXSIGS :::::::::::::::::::::::::::::::::::::
      SUBROUTINE OXSIGS(E,SIGEX,SIGEXT)
C....... Inelastic cross sections for electron impact on atomic oxygen
C....... E=electron energy, SIGEX(22)=array of partial cross sections,
C....... SIGEXT=total excitation cross section, and S
      DIMENSION SIGEX(22),SO1D(7)
      DO 9 I=1,22
 9    SIGEX(I)=0.0
      !..- CROSS SECTION FOR O(1D) - New Doering cross section from JGR
      !..- p19531, 1992. Increases production by a factor of 1.13
      DATA SO1D/0.0,0.0,15.0,30.0,44.0,54.0,38.0/
C..      IF(E.GT.6.5) SIGEX(1)=3.77E-15*(1-2.0/E)**4.25/E**1.7
C..      IF(E.LE.7.3) SIGEX(1)=SO1D(NINT(E))*1.0E-18
C
C..... Old cross section of Henry
      IF(E.GT.1.96) SIGEX(1)=4E-16*(1-1.96/E)**2/E
C........ O(1S) cross section: may be double Shyn et al. JGR 1986, 13751
      IF(E.GT.4.17) SIGEX(2)=6.54E-17*(1-SQRT(4.17/E))/E
C....... 1304, 1027 A, Zipf and Erdman JGR 1985, 11088 include cascade.
C....... Direct excitation is half for  1304 (Vaughan and Doering,
C........ JGR 1986, 13755 and 1987 in press)
      IF(E.GE.10) SIGEX(3)=67.6E-17*(E-10)/E**2
C....... 989 cross section from Doering 1987 (1/2 of Zipf)
      IF(E.GE.14) SIGEX(4)=7.0E-17*(1-14/E)/SQRT(E)
      SIGEX(5)=0.38*SIGEX(4)
C....... O(5S) 1356 A Stone And Zipf Corrected By Zipf And Erdman 1985
      !..- reparameterized 1 May 92 using FITXS.FOR (PGR)
      IF(E.GT.10.0) SIGEX(6)=4.867E-12*(1.0-9.0/E)**2.67/ E**4.0
      SIGEXT=SIGEX(1)+(SIGEX(2)+SIGEX(3)+SIGEX(4)+SIGEX(5)+SIGEX(6))
      RETURN
      END
C
C
C.................... RSPRIM.FOR ..................................
C.... This routine evaluates the ionization rates for photon impact
C.... It is based on a FLIP model routine that was modified in August 
C.... 2009 for the chemical equilibrium model by P. richards. 
      SUBROUTINE PRIMPR(IJ,Z,ZOX,ZN2,ZO2,HE,SZA,TN,F107,F107A,N4S)
      IMPLICIT NONE
      INTEGER IVERT,I,IJ,IK,IPROBS,IS,K,L,LMAX,NNI,K1
      REAL EUVION,F107,F107A,F107SV,FNFAC,FREQLY,FREQSR,O2LYXS,
     >  O2SRXS,TAUN,UVFAC,ZLAM,EUV,FLUXN,LAMAX,PEPION,PEXCIT,
     >  SIGABS,SIGION,TPOT,ZFLUX
      REAL Z,ZOX,ZN2,ZO2,HE,SZA,TN,CHI,ZZ,TNJ,TAU,FLUX,HEPLS,
     >  FBSBN,DISN,TAUGAM,FLUXG,ALTG,XNSIGF,DSPECT,GL,N4S
      REAL COLUM(3),OTHPR1,OTHPR2,OTHPR3(6)
      REAL COLUMN(3),XN(3),PROB(3,6,37),XSNPLS(37),FNITE(37),CLNITE(3)
cpgr
      REAL TPROB(3,6,37)
cpgr      

      !-- common to hold the EUV and photoelectron production rates 
      COMMON/EUVPRD/EUVION(3,12),PEXCIT(3,12),PEPION(3,12),
     >   OTHPR1(6),OTHPR2(6)
      COMMON/SIGS/ZFLUX(37),SIGABS(3,37),ZLAM(37),SIGION(3,37),
     > TPOT(3,10),NNI(3),LAMAX
      COMMON/SOL/UVFAC(59),EUV
      SAVE PROB,F107SV,TPROB    !.. Save values that are only calc once

      DATA LMAX/0/, F107SV/0.0/, IPROBS/0/
      !.. Fluxes for nighttime ion production in the 37 wavelength bins of
      !.. Torr et al GRL 1979. The fluxes are set to reproduce the production
      !.. rates in Strobel et al. PSS, p1027, 1980. Note that most bins are 
      !.. set to zero and that the Strobel production rates are scaled by 
      !.. FNFAC to stabilize the O+ solution below 200 km. Note also that
      !.. the wavelengths in FNITE go from largest (#3=HI) to smallest.
      DATA FNITE/9E5,0.0,9E5,2*0.0,9E6,13*0.0,3E5,8*0.0,3E5,8*0.0/
      DATA FNFAC/1.0/

      !.. UVFAC(58) is left over from FLIP routines for compatibility
      UVFAC(58)=-1.0 
      IF(ABS((F107-F107SV)/F107).GT.0.005) THEN
        !.. update UV flux factors
        CALL FACEUV(UVFAC,F107,F107A)
        CALL FACSR(UVFAC,F107,F107A)
        
        !.. call params to get solar flux data and cross sections
        CALL PARAMS(0,LMAX)
        F107SV=F107
      ENDIF

      !..  find probability for formation of each state  ........
      IF(IPROBS.EQ.0) THEN
        CALL PROBS(0,PROB,ZLAM,LMAX,NNI)
        DO I=1,3
        DO K=1,6
        DO L=1,37
          TPROB(I,K,L)=PROB(I,K,L)
        ENDDO
        ENDDO
        ENDDO
        IPROBS=1
       ENDIF

      !... initialization of production rates. 1.0E-15 stabilizes 
      !... e density evaluation at low altitudes in CMINOR
      DO 10 IS=1,3
      DO 10 IK=1,12
         EUVION(IS,IK)=1.0E-15
 10   CONTINUE

      DISN=0.0
      DO 687 I=1,6
      OTHPR2(I)=1.0E-15
 687  OTHPR1(I)=1.0E-15
C
C........ Nighttime He+ production is calculated and stored. Attenuated to
C........ avoid excess production at low altitudes
      OTHPR1(2)= 8E-11* EXP(-1.0E-11*ZN2) *HE
      DO 786 I=1,3
 786  COLUM(I)=1.0E+25
      TNJ=TN
      XN(1)=ZOX
      XN(2)=ZO2
      XN(3)=ZN2
      ZZ=Z*1.0E+5
      CHI=SZA
C
C*****  obtain reaction rates from subr rats to get their densities

C....... determine if sun is below the horizon ...
C---- Must now do calculation for night production - Feb 93
      ALTG=(6371.0+Z)*SIN(3.1416-CHI)-6371.0
C....      IF(CHI.GT.1.57.AND.ALTG.LT.85.) RETURN
      IF(Z.GT.1500) RETURN
C
C...... get column densities for scattered light at night  &&&&&&&&
      CALL SCOLUM(IJ,0.0E0,ZZ,TNJ,XN,CLNITE)
C
C...... evaluate the neutral column density  &&&&&&&&
      CALL SCOLUM(IJ,CHI,ZZ,TNJ,XN,COLUMN)
C........ Store the column densities for the 2-Stream program
      COLUM(1)=COLUMN(1)
      COLUM(2)=COLUMN(2)
      COLUM(3)=COLUMN(3)
C
C........ O2 dissociation by Schumann-Runge UV.
C........ OTHPR1(3)= dissociation rate. OTHPR1(5)= Energy
      CALL SCHUMN(IJ,Z,ZO2,COLUMN,OTHPR1(3),OTHPR1(5))
C
C---- Calculate hv + NO ion. freq. from Lyman-a (Brasseur & Solomon)
C---- OTHPR2(2) is photodissociation of NO in the SR bands. 
C---- A small night production from scattered light is included. FREQLY
C---- varies with solar activity using Richards et al. 1994 page 8981
C---- LY_a=2.5E11 (Lean), sigi(NO)=2.0E-18 (Brasseur & Solomon page 329)
      DATA O2LYXS,O2SRXS,FREQSR /1.0E-20,1.0E-21,5.0E-6/
       FREQLY=5.0E-7*(1+4.0E-3*(0.5*(F107+F107A)-80.0))
       OTHPR2(1)=FREQLY*(EXP(-O2LYXS*COLUMN(2))
     >    +0.001*EXP(-O2LYXS*CLNITE(2)))
       OTHPR2(2)=FREQSR*(EXP(-O2SRXS*COLUMN(2))
     >    +0.001*EXP(-O2SRXS*CLNITE(2)))
C
      !..  wavelength loop begins here  ----------
      !..  TAU, TAUN = optical depth for day, night 
      HEPLS=0.0
      DO 6 L=1,LMAX
        TAU=0.
        TAUN=0.0
        DO I=1,3
          TAUN=TAUN+SIGABS(I,L)*CLNITE(I)
          TAU=TAU+SIGABS(I,L)*COLUMN(I)
        ENDDO

        !.. evaluate nighttime flux and daytime flux
        FLUXN=FNFAC*(F107/75.0)*FNITE(L)*EXP(-TAUN)
        FLUX=ZFLUX(L)*EXP(-TAU) + FLUXN
        !..WRITE(9,'(I6,1P,22E10.2)') L,COLUMN(1),COLUMN(3),TAU,EXP(-TAU),
        !..>    FLUX, ZFLUX(L),FLUXN

        !.. he+ production. He+ X-S  = 0.25 N2  X-S. HEPRDN = nite He+
        IF(ZLAM(L).LT.500.) HEPLS=HEPLS+HE*0.25*SIGION(3,L)*FLUX

        !-- hv + N -> N+ + e. ion. freq. Richards et al. JGR 1994 page 8989
        DATA XSNPLS/6*0.0,.211,10.294,11.171,10.961,11.244,11.323,12.098
     >  ,13.265,12.423,11.951,11.212,11.798,11.758,11.778,11.772,11.503
     >  ,11.016,10.578,9.556,8.15,8.302,7.298,6.413,6.399,5.192,5.725
     >  ,4.787,3.778,2.3,.878,.286/

        OTHPR2(3)=OTHPR2(3)+XSNPLS(L)*1.0E-18*FLUX*N4S

        IF(ZLAM(L).GE.600.0) THEN
          !...... calculation of total euv absorption-ionization.....
          FBSBN=FLUX*(SIGABS(3,L)-SIGION(3,L))*XN(3)
          !.. Save energy absorbed in the photodissociative process
          OTHPR1(4)=OTHPR1(4)+1.24E+4*FBSBN/ZLAM(L)
          !.. production on atomic nitrogen by dissociation
          DISN=DISN+FBSBN
          !..      IF(J.EQ.1) WRITE(6,95) L,ZLAM(L),TAU,FLUX,FBSBN,DISN,HEPLS
          !95   FORMAT(I4,F9.1,1P,22E9.1)
          !.. take into account the large n2 absorption of lyman gamma(972.54)
          IF(NINT(ZLAM(L)).EQ.975) THEN
             TAUGAM=370E-18*COLUMN(3)
             IF(TAUGAM.GT.70.0)TAUGAM=70.0
             FLUXG=UVFAC(34) *0.82E+9 *EXP(-TAUGAM)
             DISN=DISN+FLUXG*370E-18*XN(3)
           ENDIF
         ENDIF

        !***** species loop begins here *****
        DO 304 I=1,3
          XNSIGF=XN(I)*SIGION(I,L)*FLUX
          K1=NNI(I)
 
          !.. dspect=# ions formed by w-l l by ionization of k state of species i
          DO 302 K=1,K1
            DSPECT=XNSIGF*PROB(I,K,L)
            !.. store ion production rates .....
            EUVION(I,K)=EUVION(I,K)+DSPECT

            !.. calculation of ion heating rate......
            EUVION(1,10)=EUVION(1,10)+DSPECT*TPOT(I,K)

cp           write(6,'(3I4,1P,9E10.2)') I,K,L,PROB(I,K,L),TPROB(I,K,L)
cp     >     ,XN(I),SIGION(I,L),FLUX,XNSIGF,EUVION(I,K),XSNPLS(1),FNITE(1)

 302      CONTINUE
 304    CONTINUE
 6    CONTINUE
C
      !..---   wavelength loop ends here   -----------
C
C.........Store UV disoc of N2 2 atoms produced for every dissociation
      OTHPR1(1)=2.0*DISN
C........ Transfer He+ production to storage
      OTHPR1(2)=OTHPR1(2)+HEPLS

cp      WRITE(6,*) 'After 6: EUVION(1,1), EUVION(2,1), EUVION(3,1)'
cp      write(6,'(1P,9E10.2)') EUVION(1,1), EUVION(2,1), EUVION(3,1)

 777  RETURN
      END

C:::::::::::::::::::::::::::: SCOLUM ::::::::::::::::::::::::::::::::::
C.... this routine evaluates the neutral column density for O, O2, and N2
C.... see Smith & Smith JGR 1972 p 3592
C.... chi=solar zenith angle, RE & GE radius and grav constant for earth
C.... Modified by P. Richards January 2010 to eliminate need for calling
C.... the MSIS model at grazing incidence
      SUBROUTINE SCOLUM(J,CHI,Z,TN,XN,COLUMN)
      IMPLICIT NONE
      INTEGER I,J
      REAL ZG,CHI,Z,TNJ,ALTG,GE,GR,RE,RP
      REAL SH,XP,Y,ERFY2,CHAPFN,RG,HG,XG,EM,F,G,A,B,C,D,GL
      REAL TN,XI,TINF,GTN
      REAL XN(3),COLUMN(3),SN(3),M(3),DG(9),T(2),GN(3)
      DATA A,B,C,D,F,G/1.0606963,0.55643831,1.0619896,1.7245609
     >  ,0.56498823,0.06651874/
      DATA SN/0.0,0.0,0.0/
      DATA    EM   ,   M(1) , M(2) , M(3) ,  RE   , GE
     >    / 1.662E-24 ,   16. ,  32. ,  28. ,6.357E8, 980/
      DATA T,ALTG,ERFY2/0.0,0.0,0.0D0,0.0D0/
      DATA DG/9*0.0/

      DO I=1,3
        SN(I)=0.0
        COLUMN(I)=1.E+25
      ENDDO

      TNJ=TN     !.. Avoids changing Tn at grazing incidence

      IF(CHI.LT.1.5708) GO  TO 2938      !.. is sza>90.0 degrees

      !..Grazing incidence parameters 
      ALTG=(6371.0E5+Z)*SIN(3.1416-CHI)-6371.0E5
      IF(ALTG.GE.85*1.0E5) THEN
        ZG=ALTG*1.E-5

        !.. Bates-Walker temperature
        XI=(ZG-120.0)*(6357.0+120.0)/(6357.0+ZG)
        TINF=MAX(TN,500.0)   !.. Crude TINF
        GTN=MAX(TINF-(TINF-300.0)*EXP(-0.025*XI),180.0)

        !.. Neutral densities are extrapolated from altitude to grazing
        !.. altitude. Weighted average Tn and GTn is used 
        GR=GE*(RE/(RE+Z))**2   !.. gravity
        DO I=1,3
          GN(I)=XN(I)*EXP((Z-ALTG)/
     >      ((1.38E-16*(TN+GTN*2)/3)/(EM*M(I)*GR)))
        ENDDO
        !..   WRITE(88,'(6F8.2,1P,22E10.2)') Z/1.0E5,ZG,CHI*180/3.1416,TN,GTN
        !.. >      ,TNJ,((XN(I),GN(I),SN(I)),I=1,3)
        !.. Make sure that grazing incidence density is not lower than MSIS
        !.. This is for using non MSIS densities like CTIPe
        TNJ=GTN
        DO I=1,3
          SN(I)=GN(I)
          IF(SN(I).LT.XN(I)) SN(I)=XN(I)
        ENDDO
      ELSE
        RETURN
      ENDIF
      !.. sn(1)=o , sn(2)=o2 , sn(3)=n2 , tnj=tn,  gr=gravity, rp=distance
      !.. to pt p, sh=scale height, rg=distance to pt g, hg=scale height at g

 2938 CONTINUE
      GR=GE*(RE/(RE+Z))**2
      RP=RE+Z
      !.. Calculate column densities for each species
      DO 10 I=1,3
        SH=(1.38E-16*TNJ)/(EM*M(I)*GR)
        XP=RP/SH
        Y=SQRT(0.5*XP)*ABS(COS(CHI))
        IF(Y.GT.100.0) WRITE(6,100) I,Z/1.0E5,CHI*57.3,TNJ,EM,M(I),GR,RP
 100    FORMAT('WARNING, Y IN COLUMN(I) > 100',I4,1P,9E10.2)
        IF(Y.GE.8) ERFY2=F/(G+Y)
        IF(Y.LT.8) ERFY2=(A+B*Y)/(C+D*Y+Y*Y)

 4      IF(CHI.GT.1.5708)GO  TO 2
        CHAPFN=SQRT(0.5*3.1416*XP)*ERFY2
        COLUMN(I)=XN(I)*SH*CHAPFN
        GO TO 10

 2      RG=RP*SIN(3.1416-CHI)
        HG=1.38E-16*TNJ/
     >    (EM*M(I)*980.*(6371.E5/(6371.E5+ALTG))**2)
        XG=RG/HG
        COLUMN(I)=SQRT(0.5*3.1416*XG)*HG*(2.0*SN(I)-XN(I)*ERFY2)
 10   CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PARAMS(ISW,LMAX)
C........ this program determines the cross sections, solar fluxes, and
C........ given in m. torr et al g.r.l 1979 p771, table 2 and 3. but
C........ the longer wavelengths come first
      IMPLICIT NONE
      INTEGER  I,IN,IS,ISW,J,L,LAMAX,LMAX,NNI
      REAL EUV,FFAC, SIGABS,SIGION,TPOT,UVFAC,ZFLUX,ZLAM
      REAL  X1(111),X2(111),X3(18),ZLX(37),ZFX(37)
      COMMON/SIGS/ZFLUX(37),SIGABS(3,37),ZLAM(37),SIGION(3,37),
     > TPOT(3,10),NNI(3),LAMAX
      COMMON/SOL/UVFAC(59),EUV

C
C....... ionization potentials for o,o2 and n2 see kirby et al note the
C....... o2 2(pi)u and 2(sigma-)u , and n2 f2(sigma)u pots are guesses
C....... the sixth n2 potential is for dissociation
      DATA X3/13.6,16.9,18.6,28.5,40.,0.0,12.1,16.5,18.2,20.,
     > 25.,0.,15.6,16.7,18.8,25.,29.,37./
C........ wavelength data. average is taken for bands
      DATA ZLX/1025.,1031.91,1025.72,975.,977.02,925.,875.,825.,775.,
     > 789.36,770.41,765.15,725.,703.36,675.,625.,629.73,609.76,575.,
     > 584.33,554.31,525.,475.,465.22,425.,375.,368.07,325.,303.78,
     > 303.31,275.,284.15,256.3,225.,175.,125.,75./
C........ fluxes from table 3. these are for 74113. just replace this data
C........ for other years in the table. note!!!! flux doubled for lambda<250
C........ shortest wavelenghts have been tripled
      DATA ZFX/2.4665,2.1,3.5,1.4746,4.4,3.,3.537,1.625,.758,.702,
     > .26,.17,.141,.36,.23,.342,1.59,.53,.357,1.27,.72,.452,.285,
     > .29,.383,.314,.65,.965,6.9,.8,1.679,.21,.46,3.1,4.8,.45,1.2/
C........ absorption cross sections -- o first ,o2, then n2
      DATA X1/5*0.0,1.315,4.554,3.498,5.091,3.749,3.89,4,10.736,11.46
     > ,17.245,13.365,13.4,13.4,13.024,13.09,12.59,12.059,12.127,11.93
     > ,11.496,9.687,9.84,8.693,7.7,7.68,6.461,7.08,6.05,5.202,3.732
     > ,1.839,.73,  1.346,1.0,1.63,21.108,18.73,12.817,8.562,16.631
     > ,22.145,26.668,18.91,20.8,28.535,27.44,21.919,26.017,32.06
     > ,28.07,26.61,22.79,26.04,24.606,23.101,21.91,20.31,18.118
     > ,18.32,17.438,16.81,16.8,14.387,15.79,13.37,10.9,7.509,3.806
     > ,1.316,  3*0.0,50.988,2.24,9.68,20.249,16.992,33.578,16.487
     > ,14.18,120.49,24.662,26.54,31.755,23.339,23.37,22.79,22.787
     > ,22.4,24.13,24.501,23.471,23.16,21.675,16.395,16.91,13.857
     > ,11.7,11.67,10.493,10.9,10.21,8.392,4.958,2.261,0.72/
C....... ionization cross sections 
      DATA X2/5*0.0,1.315,4.554,3.498,5.091,3.749,3.89,4,10.736,11.46
     > ,17.245,13.365,13.4,13.4,13.024,13.09,12.59,12.059,12.127,11.93
     > ,11.496,9.687,9.84,8.693,7.7,7.68,6.461,7.08,6.05,5.202,3.732
     > ,1.839,.73,  .259,0.0,1.05,13.94,15.54,9.374,5.494,6.413,10.597
     > ,10.191,8.47,11.72,23.805,23.75,21.306,24.937,31.1,26.39
     > ,26.61,22.79,26.04,24.606,23.101,21.91,20.31,18.118
     > ,18.32,17.438,16.81,16.8,14.387,15.79,13.37,10.9,7.509,3.806
     > ,1.316,  8*0.0,14.274,8.86,8.5,65.8,15.06,25.48,29.235
     > ,23.339,23.37,22.79,22.787
     > ,22.4,24.13,24.501,23.471,23.16,21.675,16.395,16.91,13.857
     > ,11.7,11.67,10.493,10.9,10.21,8.392,4.958,2.261,0.72/
C
      NNI(1)=5
      NNI(2)=5
      NNI(3)=6
      LMAX=37
      IF(ISW.NE.0) WRITE(17,95)
 95   FORMAT(/5X,'EUV fluxes, Photoabsorption, and Photoionization ',
     >  'Cross sections',
     > /4X,'I',5X,'lam',5X,'flux',4X,'sigaOX',3X,'sigaO2'
     > ,3X,'sigaN2',3X,'sigiOX',3X,'sigiO2',3X,'sigiN2',3X,'UVfac')
C
      DO 20 L=1,LMAX
      ZLAM(L)=ZLX(L)
      FFAC=UVFAC(LMAX+1-L)
      IF(ZFX(L).LT.100) ZFLUX(L)=ZFX(L)*1.0E+9*FFAC
      !..- setting up ionization potentials
      IF(L.LE.6)THEN
         TPOT(1,L)=X3(L)
         TPOT(2,L)=X3(6+L)
         TPOT(3,L)=X3(12+L)
      ENDIF
      !..- setting up cross sections
      DO 10 IS=1,3
      IN=LMAX*(IS-1)+L
      SIGABS(IS,L)=X1(IN)*1.0E-18
      SIGION(IS,L)=X2(IN)*1.0E-18
      IF(SIGABS(IS,L).LT.SIGION(IS,L)) SIGABS(IS,L)=SIGION(IS,L)
 10   CONTINUE
C
      IF(ISW.EQ.0) GO TO 20
      WRITE(17,90) L,ZLAM(L),ZFLUX(L),(SIGABS(I,L),I=1,3)
     > ,(SIGION(I,L),I=1,3),FFAC
 20   CONTINUE
C
      IF(ISW.EQ.0) RETURN
      WRITE(17,94)
 94   FORMAT(/5X,' Ionization potentials for O, O2, N2'
     > ,/2X,'4S   2D   2P   4P   2P*  -   X2   a+A  b4   B2   dis  -'
     > ,'  X2   A2   B2   C2   F2   2s')
 60   WRITE(17,91) ((TPOT(I,J),J=1,6),I=1,3)
C
      RETURN
 90   FORMAT(1X,I4,F9.2,1P,22E9.2)
 91   FORMAT(22F5.1)
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PROBS(ISW,PROB,ZLAM,LMAX,NNI)
C.... program for finding branching ratios (probabilities for various ion
C.... and molecular states) of o,o2,n2
C.... ---refs--- m. torr et al grl 1979 page 771, kirby et al atomic data
C.... and nuclear tables 1979 23,page 63
      IMPLICIT NONE
      INTEGER I,IS,ISW,J,L,LL,LMAX,NNI(3)
      REAL YO(37,5),PROB(3,6,37),ZLAM(37),SUM
C...... coefficients of o ionization cross sections from torr et al
C..... table 2
      DATA YO/.19,.486,.952,1.311,1.539,1.77,1.628,1.92,1.925,2.259
     > ,2.559,2.523,3.073,3.34,3.394,3.421,3.65,3.92,3.62,3.61,3.88,4.25
     > ,5.128,4.89,6.739,4.0,3.89,3.749,5.091,3.498,4.554,1.315,5*0.0
     > ,.206,.529,1.171,1.762,2.138,2.62,2.325,2.842,2.849,3.446,3.936
     > ,3.883,4.896,5.37,5.459,5.427,5.67,6.02,5.91,6.17,6.29,6.159
     > ,11.453,6.57,3.997,12*0.0, .134,.345,.768,1.144,1.363,1.63
     > ,1.488,1.92,1.925,2.173,2.558,2.422,2.986,3.22,3.274,3.211,3.27
     > ,3.15,3.494,3.62,3.23,2.956,0.664,14*0.0,  .062,.163,.348,.508
     > ,.598,.71,.637,.691,.693,.815,.787,.859,.541,24*0.0, .049,.13
     > ,.278,.366,.412,.35,.383,.307,.308,28*0.0/
C
C....... production of o states from torr et al table 2 (yo array)
C....... need to reverse order of yo to correspond with lambda
      DO 10 L=1,LMAX
      LL=LMAX+1-L
      SUM=YO(LL,1)+YO(LL,2)+YO(LL,3)+YO(LL,4)+YO(LL,5)
      DO 20 I=1,5
      PROB(1,I,L)=0.0
 20   IF(SUM.NE.0.0) PROB(1,I,L)=YO(LL,I)/SUM
 10    CONTINUE
C
C....... call separate subroutines for o2 and n2 probabilities
      DO 30 L=1,LMAX
      CALL PROBO2(1,L,ZLAM(L),PROB,NNI(2))
      CALL PROBN2(1,L,ZLAM(L),PROB,NNI(3))
 30   CONTINUE
C
      IF(ISW.EQ.0) RETURN
      WRITE(17,95)
 95   FORMAT(/5X,' Photoionization branching ratios for O, O2, N2'
     > ,/3X,'Lam    4S   2D   2P   4P   2P*   -   X2   a+A  b4   B2 '
     > ,'  dis   -  X2   A2   B2   C2   F2   2s')
      DO 50 L=1,LMAX
 50   WRITE(17,90) ZLAM(L),((PROB(IS,J,L),J=1,6),IS=1,3)
 90   FORMAT(F8.2,22F5.2)
       RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE PROBN2(ISW,L,ZLAM,PROB,JPTS)
C...... the n2 probabilities are taken from kirby et al tables b and c
C...... the yield of n+ is determined first then the remaining portion
C...... of the cross section is distributed amongst the n2+ ion states
C...... (x,a,b). the dissociation yield is divided between the 3 higher
C...... energy states according to wight et al. j.phys. b, 1976
C...... the 2 other states of kirby et al are not included
C
      IMPLICIT NONE
      INTEGER I,IPTS,ISW,J,JPTS,L
      REAL A(6),B(6),X(14),Y(14,6),PROB(3,6,37),SUM,YIELD,YLAM,ZLAM
      DATA IPTS/14/
      DATA X/50.,210.,240.,280.,300.,332.,428.,500.,600.,660.,660.01,
     > 720.,747.,796./
      DATA Y/5*.32,.3,.46,.404,.308,.308,.308,.42,1.,1.,5*.55,
     > .52,.46,.506,2*.589,.692,.58,2*.0,5*.13,.12,.08,
     > .09,.103,.103,4*0.0,3*0.0,.05,.1,.15,.83,1.,6*.0,3*.0,
     > .3,.4,.79,.17,7*.0,3*1.,.65,.5,.06,8*.0/
C
C...... if zlam is too big set equal to x(max)
      YLAM=ZLAM
      !.. Prevent divide by zero
      IF(ZLAM.GT.X(14)) YLAM=X(14)-1
      IF(ZLAM.LT.X(1)) YLAM=X(1)+1

       YIELD=0.0
C...... determine yield of n+, and store in prob array
       CALL YLDISS(1,YLAM,YIELD)
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(YLAM.GT.X(I).AND.YLAM.LE.X(I+1))  GO TO 20
 10   CONTINUE
 20   SUM=0.0
C...... fit straight line between points
      DO 30 J=1,JPTS
      A(J)=(Y(I+1,J)-Y(I,J))/(X(I+1)-X(I))
      B(J)=Y(I,J)-A(J)*X(I)
 30   CONTINUE
C...... determine probabilities of n2+ states
      DO 40 J=1,JPTS
      IF(J.LE.3) PROB(3,J,L)=(A(J)*YLAM+B(J))*(1-YIELD)
      IF(J.GT.3) PROB(3,J,L)=(A(J)*YLAM+B(J))*YIELD
      SUM=SUM+PROB(3,J,L)
 40   CONTINUE
C
      IF(SUM.EQ.0.0) RETURN
C....... normalise probabilities
      DO 50 J=1,JPTS
 50   PROB(3,J,L)=PROB(3,J,L)/SUM
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE YLDISS(ISW,ZLAM,YIELD)
C..... determination of dissociative yield of n+, refer to kirby et al
C..... page 66 and table b
      IMPLICIT NONE
      INTEGER ISW,I,IPTS
      REAL X(9),Y(9),ZLAM,YIELD
      DATA IPTS/9/
      DATA X/50.,210.,240.,302.,387.,477.,496.,509.,2000./
      DATA Y/.36,.36,.346,.202,.033,.041,.024,0.0,0.0/
C
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(ZLAM.GE.X(I).AND.ZLAM.LT.X(I+1))  GO TO 20
 10   CONTINUE
 20   IF(ZLAM.GT.387.AND.ZLAM.LT.477) GO TO 40
C....... linear interpolation
      YIELD=(ZLAM-X(I))/(X(I+1)-X(I))*(Y(I+1)-Y(I))+Y(I)
      GO TO 30
C...... parabolic interpolation see formula page 66 kirby et al
 40   YIELD=.0329+8.13E-6*(ZLAM-442)**2
 30   RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE PROBO2(ISW,L,ZLAM,PROB,JPTS)
C....... o2 branching ratios are taken from kirby et al table d
C....... columns 4 & 9 are combined. columns 5,6,7&8 are combined
      IMPLICIT NONE
      INTEGER ISW,I,IPTS,J,JPTS,L
      REAL A(5),B(5),X(20),Y(20,5),PROB(3,6,37),SUM,YLAM,ZLAM
      DATA IPTS,X/20,304.,323.,454.,461.,504.,537.,556.,573.,584.,598.
     >  ,610.,637.,645.,662.,684.,704.,720.,737.,774.,1026./
      DATA Y/.365,.374,.432,.435,.384,.345,.356,.365,.306,.23,.235,.245,
     > .34,.27,.482,.675,.565,.565,1.,1.,.205,.21,.243,.245,.27,.29,
     > .23,.27,.33,.295,.385,.35,.305,.385,.518,.325,.435,.435,2*.0,
     > .125,.124,.12,.12,.126,.13,.225,.216,.21,.375,.305,.37,.33,.345,
     > 6*.0,.055,.167,.11,.105,.194,.234,.189,.149,.155,.103,.075,.036
     > ,.025,7*0.,.25,.125,.095,.95,.026,15*0./
C
C...... if zlam is too big set equal to x(max)
      YLAM=ZLAM
C...... if zlam is outside range of data values set equal to max or min
      IF(ZLAM.GT.X(20)) YLAM=X(20)
      IF(ZLAM.LE.X(1)) YLAM=X(1)+1.E-3
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(YLAM.GT.X(I).AND.YLAM.LE.X(I+1))  GO TO 20
 10   CONTINUE
 20   SUM=0.0
C
      DO 30 J=1,JPTS
      A(J)=(Y(I+1,J)-Y(I,J))/(X(I+1)-X(I))
      B(J)=Y(I,J)-A(J)*X(I)
      SUM=SUM+A(J)*YLAM+B(J)
 30   CONTINUE
C
      DO 40 J=1,JPTS
      PROB(2,J,L)=(A(J)*YLAM+B(J))/SUM
 40   CONTINUE
C
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SCHUMN(J,Z,ZO2,COLUMN,SCHUPR,SCHUHT)
C......... production of o(1d) by schumann-runge bands
C......... The fluxes are from Torr et al. GRL 1980 p6063. Scaling is
C......... done using UVFAC which may be set according to F10.7 cm flux
C......... may be done in FACEUV
      IMPLICIT NONE
      INTEGER J,JTI,LMAX,LSR
      REAL Z,ZO2,SCHUPR,SCHUHT,HSRX,FLD,EUV,SRXSCT,UVFAC
      COMMON/SOL/UVFAC(59),EUV
      REAL COLUMN(3),SRFLUX(8),SRXS(8),SRLAM(8)
      DATA SRFLUX/2.4,1.4,.63,.44,.33,.17,.12,.053/
      DATA SRXS/.5,1.5,3.4,6,10,13,15,12/
      DATA SRLAM/1725,1675,1625,1575,1525,1475,1425,1375/
C
C........ lmax=# of lambdas in sub. primpr: schuht=heating: schupr=o(1d) prod
      LMAX=37
C
      DO 505 LSR=1,8
C......... photoabsorption cross section
      SRXSCT=1.0E-18*SRXS(LSR)
      HSRX=SRXSCT*COLUMN(2)
      IF(HSRX.GT.70)HSRX=70
C........ attentuated solar flux
      FLD=UVFAC(LMAX+LSR)*1.E+11*SRFLUX(LSR)*EXP(-HSRX)
C............ neutral heating SCHUHT and photodissociation rate SCHUPR
      SCHUHT=SCHUHT+1.24E+4*(FLD*SRXSCT)*ZO2/SRLAM(LSR)
      SCHUPR=SCHUPR+FLD*SRXSCT
      !IF(JTI.EQ.0) WRITE(24,90) LSR,SRXSCT,FLD,SCHUPR,COLUMN(2),FLD,
      !>   UVFAC(LMAX+LSR)
 505  CONTINUE
C
      SCHUPR=ZO2*SCHUPR
 90   FORMAT(2X,I5,1P,9E9.1)
      JTI=1
      RETURN
      END
C:::::::::::::::::::::::::::::::: FACEUV :::::::::::::::::::::::
      SUBROUTINE FACEUV(UVFAC,F107,F107A)
C----- This routine uses the EUV scaling from Richards et al.[1994]
C----- The EUVAC flux model is based on the F74113 solar reference
C----- spectrum and Hinteregger's scaling factors. This subroutine
C----- just provides the scaling factors as a function of the proxy
C----- (F107+F107A)/2
      IMPLICIT NONE
      INTEGER I
      REAL UVFAC(59),HFG200(37),A,B,C,F107,F107A,F107AV
      DATA HFG200/2.202,1.855,2.605,3.334,1.333,17.522,4.176,4.0
     >  ,1.4,3.694,1.791,5.385,1.889,1.899,3.427,2.051,1.392,1.619
     >  ,1.439,2.941,1.399,2.416,1.512,1.365,1.570,1.462,2.537,1.393
     >  ,1.572,1.578,1.681,1.598,1.473,1.530,1.622,1.634,1.525/

      !--  Test to see if need to scale - see DATRD2 subroutine      
      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3) THEN
         !........... EUV scaling
         F107AV=(F107+F107A)*0.5
         DO 50 I=1,37
         A=(HFG200(I)-1)/120.0
         B=1-A*80.0
         UVFAC(I)=A*F107AV+B
         IF(UVFAC(I).LT.0.8) UVFAC(I)=0.8
 50      CONTINUE
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE FACSR(UVFAC,F107,F107A)
C........ The Schumann-Runge factors are scaled according to F10.7
C........ from Torr et al. GRL 1980 p6063
      IMPLICIT NONE
      INTEGER I,LSR
      REAL UVFAC(59),SRFLUX(8),SRA(8),SRB(8),F107,F107A
      !............. Schumann-Runge scaling
      DATA SRFLUX/2.4,1.4,.63,.44,.33,.17,.12,.053/
      !...... first two SRA and SRB values out of order in Marsha's paper
      DATA SRA/25.5,20.7,13.2,11.6,11.3,7.86,7.68,4.56/
      DATA SRB/222.,129.,53.4,36.0,25.0,11.3,6.35,2.05/
C
C----  Test to see if need to scale - see DATRD2 subroutine      
      !IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3) THEN
C
         DO 505 I=38,50
            LSR=I-37
            UVFAC(I)=1.0
            IF(LSR.LE.8)
     >      UVFAC(I)=(SRA(LSR)*1.0E7*F107+SRB(LSR)*1.0E9)/SRFLUX(LSR)
     >           /1.0E11
 505     CONTINUE
      !ENDIF
      RETURN
      END
