#include "dims.h"
      SUBROUTINE QINITE
      use cons_module,only: len1,len2,len3,kmax,kmaxm1,kmaxp1,expz,
     |  rmassinv,p0,boltz,expzmid_inv,expzmid
      implicit none
C     ****
C     ****      CALCULATE BACKGROUND IONIZATION RATES
C     ****
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "crates_tdep.h"
#include "cmpbnd.h"
#include "mwt.h"
#include "phys.h"
!
! Local:
      real :: SA(3,3),SI(3,3),AL(3)
      integer :: i,n,npsk,k,nps1k,nps2k,ntk,nmsk,m,nnvk,nqo2pk,nqopk,
     |  nqn2pk,nqnpk,nqnopk,nnvo2k,npnok
!
      DATA
     1SA/1.6E-18,0.,0.,22.E-18,10.24E-18,23.11E-18,16.E-18,8.4E-18,
     211.61E-18/,
     3SI/1.0E-18,0.,0.,22.E-18,10.24E-18,23.11E-18,16.E-18,8.4E-18,
     411.61E-18/,
C     ****
C     ****    MOD FOR QINITE (COMPUTES BACKGROUND IONIZATION RATES)
C     ****
C    5AL/5.E8,5.E7,5.E7/
C    5AL/5.E7,5.E6,5.E6/
!    5AL/2.0E7,2.0E6,2.0E6/  ! snoe mod
     5AL/5.E4,5.E3,5.E3/     ! tgcm13mt
C    5AL/1.0E7,1.0E6,1.0E6/
C    5AL/1.5E7,1.5E6,1.5E6/
C    5AL/8.0E6,8.0E5,8.0E5/
C    5AL/5.5E6,5.5E5,5.5E5/
C    5AL/5.E6,5.E5,5.E5/
C    5AL/5.E5,5.E4,5.E4/
C    5AL/1.E5,1.E4,1.E4/
C    5AL/5.E4,5.E3,5.E3/
C     ****
C     ****      S1 = PSI1, S2 = PSI2, S3 = PSI3          (K)
C     ****
      DO 1 N=1,2
      NPSK= NJ+NPS+(N-1)*KMAXP1
C     ****      LEVELS 2 THRU KMAXP1
      K=(2-N)*KMAXP1+2
      DO 1 I=1,LEN2
      S2(I,K)=.5*(F(I,NPSK)+F(I,NPSK+1))
    1 CONTINUE
C     ****      LEVEL 1
      NPS1K=NJ+NPS
      NPS2K=NPS1K+KMAXP1
      DO 2 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))
    2 CONTINUE
C     ****      S3 = PSI3               (K)
      DO 3 I=1,LEN3
      S3(I,1)=1.-S1(I,1)-S2(I,1)
    3 CONTINUE
C     ****
C     ****      S4 =T                   (K)
C     ****
      NTK=NJ+NT
C     ****      LEVELS 2 THRU KMAX
      DO 4 I=LEN1+1,LEN2
      S4(I,1)=.5*(F(I,NTK)+F(I,NTK-1))
    4 CONTINUE
C     ****      LEVELS 1 AND KMAXP1
      DO 5 I=1,LEN1
      S4(I,1)=F(I,NTK+KMAX)
      S4(I,KMAXP1)=F(I,NTK+KMAXM1)
    5 CONTINUE
C     ****
C     ****      S5 =N*MBAR              (K)
C     ****
C     ****      S5 = expz               (K)
      DO 6 K=1,KMAX
      DO 6 I=1,LEN1
      S5(I,K)=expzmid_inv*expz(K)
    6 CONTINUE
      DO 7 I=1,LEN1
      S5(I,KMAXP1)=expzmid*expz(KMAX)
    7 CONTINUE
C     ****      S3 = N*MBAR
      NMSK=NJ+NMS
      DO 8 I=1,LEN3
      S5(I,1)=p0*S5(I,1)*F(I,NMSK)/(boltz*S4(I,1))
    8 CONTINUE
C     ****
C     ****      CALCULATE IONIZATION OF O2, O, N2 IN S6, S7, S8
C     ****
C     ****      INITIALIZATION  S6 = S7 = S8 = 0.
      DO 9 I=1,3*LEN3
      S8(I,1)=0.
    9 CONTINUE
C     ****      SUMMATION OVER WAVE LENGTH
      DO 10 N=1,3
C     ****      S9 = TAU
      DO 11 I=1,LEN3
      S9(I,1)=0.
   11 CONTINUE
C     ****      SUMMATION OVER SPECIES
      DO 12 M=1,3
      NNVK=NNVO2+(M-1)*KMAXP1
      DO 12 I=1,LEN3
      S9(I,1)=S9(I,1)+SA(M,N)*F(I,NNVK)
   12 CONTINUE
C     ****      S9 = EXP(-TAU) = EXP(-S9)
      DO 13 I=1,LEN3
      S9(I,1)=EXP(-S9(I,1))
   13 CONTINUE
C     ****      ADD O2, O, N2 IONIZATION TO S6, S7, S8
      DO 14 M=1,3
      K=(3-M)*KMAXP1+1
      DO 14 I=1,LEN3
      S8(I,K)=S8(I,K)+AL(N)*SI(M,N)*S3(I,K)*S9(I,1)*rmassinv(M)
   14 CONTINUE
   10 CONTINUE

!     call addfsech('TAU' ,' ',' ',s9 ,zimxp,zkmxp,zkmxp,j)
!     call addfsech('QBO2',' ',' ',s6 ,zimxp,zkmxp,zkmxp,j)
!     call addfsech('QBO1',' ',' ',s7 ,zimxp,zkmxp,zkmxp,j)
!     call addfsech('QBN2',' ',' ',s8 ,zimxp,zkmxp,zkmxp,j)

C     ****
C     ****      CALCULATE CONTRIBUTIONS TO NQO2P, NQOP, NQN2P, NQNP
C     ****
      NQO2PK=NQO2P
      NQOPK=NQOP
      NQN2PK=NQN2P
      NQNPK=NQNP
!
! 4/28/00: Original snoe mods had NQO2P+1.E+2. Also, NQOP, NQN2P, and 
!          NQNP were missing. The addition of 100 to NQO2P was resulting
!          in incorrect O2+ above about zp=0 (elden.f). This also
!          affected a TE problem, which was later solved with use of
!          DIPMIN in lamdas and settei.
!
      DO 15 I=1,LEN3
        F(I,NQO2PK)=F(I,NQO2PK)+0.67*S6(I,1)*S5(I,1)
        F(I,NQOPK) =F(I,NQOPK) +(0.33*S6(I,1)+S7(I,1))*S5(I,1)
        F(I,NQN2PK)=F(I,NQN2PK)+0.86*S8(I,1)*S5(I,1)
        F(I,NQNPK) =F(I,NQNPK) +0.14*S8(I,1)*S5(I,1)
   15 CONTINUE
C     ****
C     ****      CALCULATE NO IONIZATION AND ADD TO NQNOP
C     ****
      NQNOPK=NQNOP
      NNVO2K=NNVO2
      NPNOK=NJ+NPNO
C     ****      LEVEL 1
      DO 16 I=1,LEN1
      F(I,NQNOPK)=F(I,NQNOPK)+BETA9N(I,1)*F(I,NPNOK)
     1*S5(I,1)/RMNO
   16 CONTINUE
C     ****      LEVELS 2 THRU KMAXP1
      DO 17 I=LEN1+1,LEN3
      F(I,NQNOPK)=F(I,NQNOPK)+BETA9N(I,1)*.5*(F(I,NPNOK)+
     1F(I,NPNOK-1))*S5(I,1)/RMNO
   17 CONTINUE

!     call addfsech('QO2P',' ',' ',f(1,nqo2p),zimxp,zkmxp,zkmxp,j)
!     call addfsech('QOP' ,' ',' ',f(1,nqop ),zimxp,zkmxp,zkmxp,j)
!     call addfsech('QN2P',' ',' ',f(1,nqn2p),zimxp,zkmxp,zkmxp,j)
!     call addfsech('QNP' ,' ',' ',f(1,nqnp ),zimxp,zkmxp,zkmxp,j)
!     call addfsech('QNOP',' ',' ',f(1,nqnop),zimxp,zkmxp,zkmxp,j)

      RETURN
      END
C