!
      module qjion_module
!
! Heating due to ion chemistry.
!
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine qjion(tn,o2,ox,n2,o2p,op,n2p,nplus,nop,xiop2p,xiop2d,
     |  n4s,n2d,no,ne,o1,barm,cp,qic,qphoto,lev0,lev1,lon0,lon1,lat)
!
      use qrj_module,only:
     |  qo2p,   ! o2+ ionization            ! F(NQO2P)
     |  qop,    ! o+  ionization            ! F(NQOP)
     |  qn2p,   ! n2+ ionization            ! F(NQN2P)
     |  qnop,   ! no+ ionization            ! F(NQNOP)
     |  qnp,    ! n+  ionization            ! F(NQNP)
     |  qtotal  ! total heating             ! F(NQ)
      use cons_module,only: 
     |  avo,p0,expz,expzmid,expzmid_inv,boltz,evergs,
     |  rmi_o2=>rmassinv_o2  ,rmi_n2=>rmassinv_n2,rmi_n2d=>rmassinv_n2d,
     |  rmi_n4s=>rmassinv_n4s,rmi_no=>rmassinv_no,rmi_o1=>rmassinv_o1
      use chemrates_module,only:
     |  rk1,rk2,rk3,rk4,rk5,rk6,rk7,rk8,rk9,rk10,rk16,rk17,rk18,rk19,
     |  rk20,rk22,rk23,rk24,rk25,rk26,rk27,rk28,ra1,ra2,ra3
      use addfld_module,only: addfld
      use diags_module,only: mkdiag_QIC, mkdiag_QPHOTO
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,          ! neutral temperature    (f(:,nj+nt  ))
     |  o2,          ! molecular oxygen       (f(:,nj+nps ))
     |  ox,          ! ox (mmr)               (f(:,nj+nps2))
     |  n2,          ! n2 (mmr)               (f(:,nj+npn2))
     |  o2p,         ! o2+                    (f(:,nj+no2p))
     |  op,          ! o+                     (f(:,nj+nop ))
     |  n2p,         ! n2+                    (f(:,nj+nn2p))
     |  nplus,       ! nplus                  (f(:,nj+nnp ))
     |  nop,         ! no+                    (f(:,nj+nnop))
     |  xiop2p,      ! o+(2p)                 xiop2p
     |  xiop2d,      ! o+(2d)                 xiop2d
     |  n4s,         ! n(4s) (mmr)            (f(:,nj+npn4s))
     |  n2d,         ! n(2d) (mmr)            (f(:,nj+npn2d))
     |  no,          ! no (mmr)               (f(:,nj+npno))
     |  ne,          ! electron density (cm3) (f(:,nj+ne  ))
     |  o1,          ! atomic oxygen (mmr)    (f(:,nj+npo1))
     |  barm,        ! mean molecular weight  (f(:,nj+nms ))
     |  cp           ! mean molecular weight  (f(:,nj+nms ))
!
! Output args:
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  qphoto,      ! photo electron heating of neutral gas ! s2
     |  qic          ! ion chemistry heating of neutral gas  ! s3
!
! Local:
      integer :: i,k,i0,i1,nk,nkm1
      real,parameter :: aureff = 0.05
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  xnmbari,     ! p0e(-z)*barm/(boltz*T) at interfaces  ! s11
     |  xnmbarm,     ! p0e(-z)*barm/(boltz*T) at midpoints   ! s12
     |  qt,          ! total ionization rate over 5 species  ! s1
     |  diag_qphoto, ! QPHOTO diagnostic (deg K/day)
     |  diag_qic     ! QIC diagnostic (deg K/day)
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! p0e(-z)*mbar/kT at midpoints and interfaces (boundaries are handled
! slightly differently than xnmbar variables in the 3-d vars in fields.F,
! so find them locally here):
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          xnmbarm(k,i) = p0*expz(k)*.5*(barm(k,i)+barm(k+1,i))/
     |      (boltz*tn(k,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          xnmbari(k,i) = p0*expzmid_inv*expz(k)*barm(k,i)/
     |      (boltz*.5*(tn(k,i)+tn(k-1,i)))
        enddo ! k=lev0,lev1-1
        xnmbari(lev0,i) = p0*expzmid_inv*expz(lev0)*barm(lev0,i)/ ! bottom
     |    (boltz*.5*(3.*tn(lev0,i)-tn(lev0+1,i)))
        xnmbari(lev1,i) = p0*expzmid*expz(lev1-1)*barm(lev1,i)/     ! top
     |    (boltz*.5*(3.*tn(lev1-1,i)-tn(lev1-2,i)))
      enddo ! i=lon0,lon1

!     call addfld('XNMBARM',' ',' ',xnmbarm(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('XNMBARI',' ',' ',xnmbari(:,lon0:lon1),
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!
! Total ionization rate over 5 species:
!
      qt = 0. ! whole array op
      do i=lon0,lon1
        do k=lev0,lev1
          qt(k,i) = qt(k,i)+qo2p(k,i,lat)+qop(k,i,lat)+qn2p(k,i,lat)+ ! s1
     |                      qnop(k,i,lat)+qnp(k,i,lat)
          qphoto(k,i) = qt(k,i)*aureff*35.*avo*evergs/xnmbari(k,i)    ! s2
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
! Ion chemstry heating:
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          qic(k,i) = (avo*(o2(k,i)*rmi_o2*(rk1(k,i,lat)*op(k,i)*1.555+ ! s3
     |      (rk6*2.486+rk7*6.699)*nplus(k,i)+rk9*n2p(k,i)*3.52)+
     |      op(k,i)*(rk2(k,i,lat)*n2(k,i)*rmi_n2*1.0888+
     |      rk10*n2d(k,i)*rmi_n2d*1.45)+
     |      ox(k,i)*rmi_o1*(rk3(k,i,lat)*n2p(k,i)*0.70+
     |      rk8*nplus(k,i)*0.98)+o2p(k,i)*(rk4*n4s(k,i)*rmi_n4s*4.21+
     |      rk5*no(k,i)*rmi_no*2.813))+.5*(ne(k,i)+ne(k+1,i))*
     |      (ra1(k,i,lat)*nop(k,i)*0.854+ra2(k,i,lat)*o2p(k,i)*5.2755+
     |       ra3(k,i,lat)*n2p(k,i)*3.678)/xnmbarm(k,i))*evergs+
     |      (avo*(((rk16*3.02+rk17*0.7)*n2(k,i)*rmi_n2+
     |      rk18*o1(k,i)*rmi_o1*5.0+rk23*o2(k,i)*rmi_o2*6.56)*
     |      xiop2p(k,i)+(rk24*n2(k,i)*rmi_n2*1.33+
     |      rk25*ox(k,i)*rmi_o1*3.31+rk27*4.87*o2(k,i)*rmi_o2)*
     |      xiop2d(k,i))+(0.5*(ne(k,i)+ne(k+1,i))*((rk19(k,i,lat)*5.0+
     |      rk20(k,i,lat)*1.69)*xiop2p(k,i)+
     |      rk26(k,i,lat)*3.31 *xiop2d(k,i))-rk22*1.69*xiop2p(k,i)-
     |      rk28*3.33*xiop2d(k,i))/xnmbarm(k,i))*evergs

          if (qic(k,i) < 1.e-30) qic(k,i) = 1.e-30

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Add to qtotal:
!

      do i=lon0,lon1
        do k=lev0,lev1-2
          qtotal(k+1,i,lat) = qtotal(k+1,i,lat)+(qphoto(k+1,i)+
     |      sqrt(qic(k,i)*qic(k+1,i)))
        enddo ! k=lev0,lev1-2
        qtotal(lev0,i,lat) = qtotal(lev0,i,lat)+(qphoto(lev0,i)+  ! bottom
     |    sqrt(qic(lev0,i)**3/qic(lev0+1,i)))
        qtotal(lev1,i,lat) = qtotal(lev1,i,lat)+(qphoto(lev1,i)+  ! top
     |    sqrt(qic(lev1-1,i)**3/qic(lev1-2,i)))
      enddo ! i=lon0,lon1
!
! Output diagnostics in deg K/day:
!   QPHOTO: photo electron heating of neutral gas
!   QIC: ion chemistry heating of neutral gas
!
      do i=lon0,lon1
        do k=lev0,lev1
          diag_qphoto(k,i) = qphoto(k,i)*86400./cp(k,i)
          diag_qic   (k,i) = qic   (k,i)*86400./cp(k,i)
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
      
      call mkdiag_QPHOTO('QPHOTO',diag_qphoto,lev0,lev1,lon0,lon1,lat)
      call mkdiag_QIC   ('QIC'   ,diag_qic   ,lev0,lev1,lon0,lon1,lat)

      end subroutine qjion
!-----------------------------------------------------------------------
      end module qjion_module
