#include "dims.h"
      SUBROUTINE QJION
      use cons_module,only: len1,len2,len3,kmax,kmaxp1,kmaxm1,expz,
     |  p0,boltz,expzmid,expzmid_inv,avo,rmassinv_o2,rmassinv_o,
     |  rmassinv_n2
      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_const.h"
#include "crates_tdep.h"
#include "mwt.h"
C
      real :: AUREFF(ZKMXP)
      integer :: i,k,n,ntk,nmsk,nps1k,nps2k,no2pk,nopk,nn2pk,nnpk,nnopk,
     |  npn4sk,npn2dk,npnok,nek,npo1k,nqk,nqpk
      AUREFF=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
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,5
      NQPK=NQO2P+(N-1)*KMAXP1
      DO 4 I=1,LEN3
      S1(I,1)=S1(I,1)+F(I,NQPK)
    4 CONTINUE
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
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
      NPO1K=NJ+NPO1
      DO 6 I=1,LEN2
      S3(I,1)=(avo*(F(I,NPS1K)*rmassinv_o2*(RK1(I,1)*F(I,NOPK)*1.555+
     1(RK6*2.486+RK7*6.699)*F(I,NNPK)+RK9*F(I,NN2PK)*
     23.52)+F(I,NOPK)*(RK2(I,1)*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2*
     31.0888+RK10*F(I,NPN2DK)/RMN2D*1.45)+F(I,NPS2K)*rmassinv_o*
     4(RK3(I,1)*F(I,NN2PK)*0.70+RK8*F(I,NNPK)*0.98)+F(I,NO2PK)*
     5(RK4*F(I,NPN4SK)/RMN4S*4.21+RK5*F(I,NPNOK)/RMNO*2.813))
     6+.5*(F(I,NEK)+F(I,NEK+1))*(RA1(I,1)*F(I,NNOPK)*0.854+RA2(I,1)*
     7F(I,NO2PK)*5.2755+RA3(I,1)*F(I,NN2PK)*3.678)/S12(I,1))*1.602E-12
     1  +(avo*(((RK16*3.02+RK17*0.7)*(1.-F(I,NPS1K)-
     2  F(I,NPS2K))*rmassinv_n2+RK18*F(I,NPO1K)*rmassinv_o*5.0+
     3  RK23*F(I,NPS1K)*rmassinv_o2*6.56)*
     4  XIOP2P(I,1)+(RK24*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2*1.33
     5  +RK25*F(I,NPS2K)*rmassinv_o*3.31+RK27*4.87*F(I,NPS1K)*
     6  rmassinv_o2)*XIOP2D(I,1))+(0.5*(F(I,NEK)+F(I,NEK+1))*((RK19(I,1)
     7  *5.0+RK20(I,1)*1.69)*XIOP2P(I,1)+RK26(I,1)*3.31*XIOP2D(I,1))
     8  -RK22*1.69*XIOP2P(I,1)-RK28*3.33
     9  *XIOP2D(I,1))/S12(I,1))*1.602E-12
!     S3(I,1) = merge(S3(I,1),1.E-30,S3(I,1)-1.E-30>=0.)
      if (s3(i,1) < 1.e-30) s3(i,1) = 1.e-30
    6 CONTINUE
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
      RETURN
      END