#include "dims.h"
      SUBROUTINE QJION
      use cons_module,only: len1,len2,len3,kmax,kmaxm1,kmaxp1,
     |  rmassinv,expz,expzmid,expzmid_inv,avo,p0,boltz
      use crates_module,only: rk4,rk5,rk6,rk7,rk8,rk9,rk10,rk16,rk17,
     |  rk18,rk21,rk22,rk23,rk24,rk26,rk27
      implicit none
C     ****
C     ****     CALCULATES ION CHEMISTRY CONTRIBUTION TO NEUTRAL GAS
C     ****     HEATING AND O2 DISSOCIATION
C     ****
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
#include "crates_tdep.h"
#include "mwt.h"
!
! Local:
      real :: AUREFF(ZKMXP)
      integer :: i,k,ntk,nmsk,n,nqpk,nps1k,nps2k,no2pk,nopk,nn2pk,nnpk,
     |  nnopk,npn4sk,npn2dk,npnok,nek,nqk
!
      DATA AUREFF/ZKMXP*0.05/
C     ****
C     ****     S12 = N*MBAR =PO*EXP(-Z)*MBAR/(K*T)      (K+1/2)
C     ****
      NTK=NJ+NT-1
      NMSK=NJ+NMS-1
      DO 2 K=1,KMAX
      NTK=NTK+1
      NMSK=NMSK+1
      DO 2 I=1,LEN1
      S12(I,K)=p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz*F(I,NTK))
    2 CONTINUE
C     ****
C     ****     S11 = N*MBAR =PO*EXP(-Z)*MBAR/(K*T)      (K)
C     ****
      NTK=NJ+NT
      NMSK=NJ+NMS
C     ****     LEVELS 2 THRU KMAX
      DO 12 K=2,KMAX
      NTK=NTK+1
      NMSK=NMSK+1
      DO 12 I=1,LEN1
      S11(I,K)=p0*expzmid_inv*expz(K)*F(I,NMSK)/(boltz*.5*(F(I,NTK)+
     1F(I,NTK-1)))
   12 CONTINUE
C     ****     LEVELS 1 AND KMAXP1
      NTK=NJ+NT
      NMSK=NJ+NMS
      DO 13 I=1,LEN1
      S11(I,1)=p0*expzmid_inv*expz(1)*F(I,NMSK)/(boltz*.5*(3.*F(I,NTK)-
     1F(I,NTK+1)))
      S11(I,KMAXP1)=p0*expzmid*expz(KMAX)*F(I,NMSK+KMAX)/(boltz*.5*
     1(3.*F(I,NTK+KMAXM1)-F(I,NTK+KMAXM1-1)))
   13 CONTINUE

!     call addfsech('XNMBAR ',' ',' ',s11,zimxp,zkmxp,zkmx,j)
!     call addfsech('XNMBAR ',' ',' ',s12,zimxp,zkmxp,zkmx,j)
!
!     call addfsech('XNMBARM',' ',' ',s11,zimxp,zkmxp,zkmx,j)
!     call addfsech('XNMBARM',' ',' ',s12,zimxp,zkmxp,zkmx,j)
!
!     call addfsech('XNMBARI',' ',' ',s11,zimxp,zkmxp,zkmx,j)
!     call addfsech('XNMBARI',' ',' ',s12,zimxp,zkmxp,zkmx,j)

!     call addfsech('XNMBARI',' ',' ',s11,zimxp,zkmxp,zkmxp,j)
!     call addfsech('XNMBARM',' ',' ',s12,zimxp,zkmxp,zkmx ,j)

C     ****
C     ****     S1 = QT   (TOTAL IONIZATION RATE)      (K)
C     ****
C     ****     S1 = 0.
      DO 3 I=1,LEN3
      S1(I,1)=0.
    3 CONTINUE
C     ****     SUMMATION OVER 5 ION SPECIES
      DO 4 N = 1,7
      NQPK=NQO2P+(N-1)*KMAXP1
      DO 4 I=1,LEN3
      S1(I,1)=S1(I,1)+F(I,NQPK)
    4 CONTINUE

!     call addfsech('QTOT',' ',' ',s1,zimxp,zkmxp,zkmx,j)

C     ****
C     ****     S2 = QPHOTO     (PHOTO ELECTRON HEATING OF NEUTRAL GAS) (
C     ****
      DO 5 K=1,KMAXP1
      DO 5 I=1,LEN1
      S2(I,K)=S1(I,K)*AUREFF(K)*35.*avo*1.602E-12/S11(I,K)
    5 CONTINUE

!     call addfsech('QPHOTO',' ',' ',s2,zimxp,zkmxp,zkmxp,j)

C     ****
C     ****     S3 = QIC  (ION CHEMISTRY HEATING OF NEUTRAL GAS) (K=1/2)
C     ****
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NO2PK=NJ+NO2P
      NOPK=NJ+NOP
      NN2PK=NN2P
      NNPK=NNP
      NNOPK=NNOP
      NPN4SK=NJ+NPN4S
      NPN2DK=NJ+NPN2D
      NPNOK=NJ+NPNO
      NEK=NJ+NE
      DO 6 I=1,LEN2
      S3(I,1)=(avo*(F(I,NPS1K)*rmassinv(1)*(RK1(I,1)*F(I,NOPK)*
     |  1.555+(RK6*2.486+RK7*6.699)*F(I,NNPK)+RK9*
     |  F(I,NN2PK)*3.52)+F(I,NOPK)*(RK2(I,1)*(1.-F(I,NPS1K)-F(I,NPS2K))*
     |  rmassinv(3)*1.0888+RK10*F(I,NPN2DK)/RMN2D*1.45)+F(I,NPS2K)*
     |  rmassinv(2)*(RK3(I,1)*F(I,NN2PK)*0.70+RK8*F(I,NNPK)*0.98)+
     |  F(I,NO2PK)*(RK4*F(I,NPN4SK)/RMN4S*4.21+RK5*F(I,NPNOK)/
     |  RMNO*2.813))+.5*(F(I,NEK)+F(I,NEK+1))*(RA1(I,1)*F(I,NOPK)*0.854+
     |  RA2(I,1)*F(I,NO2PK)*5.2755+RA3(I,1)*F(I,NN2PK)*3.678)/S12(I,1))*
     |  1.602E-12+(avo*(((RK16*3.02+RK17*0.7)*(1.-F(I,NPS1K)
     |  -F(I,NPS2K))*rmassinv(3)+RK18*F(I,NPS2K)*rmassinv(2)*5.0)*
     |  XIOP2P(I,1)+(RK23*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)*
     |  1.33+RK24*F(I,NPS2K)*rmassinv(2)*3.31+RK26*4.87*
     |  F(I,NPS1K)*rmassinv(1))*XIOP2D(I,1))+(0.5*(F(I,NEK)+F(I,NEK+1))*
     |  ((RK19(I,1)*5.0+RK20(I,1)*1.69)*XIOP2P(I,1)+RK25(I,1)*3.31*
     |  XIOP2D(I,1))-(RK21*5.02+RK22*1.69)*XIOP2P(I,1)-
     |  RK27*3.33*XIOP2D(I,1))/S12(I,1))*1.602E-12
      if (s3(i,1) < 1.e-30) s3(i,1) = 1.e-30
    6 CONTINUE

!     call addfsech('QIC',' ',' ',s3,zimxp,zkmxp,zkmx,j)

C     ****
C     ****     ADD QPHOTO AND QIC TO Q    (K)
C     ****
      NQK=NQ
C     ****     LEVELS 2 THRU KMAX
      DO 7 I=1,LEN2-LEN1
      S6(I,2)=S2(I,2)+SQRT(S3(I,1)*S3(I,2))
      F(I,NQK+1)=F(I,NQK+1)+S6(I,2)
    7 CONTINUE
C     ****     LEVELS 1 AND KMAXP1
      DO 8 I=1,LEN1
      S6(I,1)=S2(I,1)+SQRT(S3(I,1)**3/S3(I,2))
      F(I,NQK)=F(I,NQK)+S6(I,1)
      S6(I,KMAXP1)=S2(I,KMAXP1)+SQRT(S3(I,KMAX)**3/S3(I,KMAXM1))
      F(I,NQK+KMAX)=F(I,NQK+KMAX)+S6(I,KMAXP1)
    8 CONTINUE

!     call addfsech('QTOTAL',' ',' ',f(1,nq),zimxp,zkmxp,zkmx,j)

      RETURN
      END
C