!
      subroutine qjnno(tn,o2,o1,o3,o1d,ne,no,n4s,n2d,barm,xnmbar,cp,
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculates additions to neutral gas heating and O2 dissociation
! due to N, NO chemistry.
!
      use cons_module,only: evergs,avo,expz,
     |  rmi_n2 =>rmassinv_n2 ,rmi_o2 =>rmassinv_o2 ,rmi_no=>rmassinv_no,
     |  rmi_n4s=>rmassinv_n4s,rmi_n2d=>rmassinv_n2d,rmi_o3=>rmassinv_o3,
     |  rmi_o1 =>rmassinv_o1
      use chemrates_module,only: beta1,beta2,beta3,beta4,beta5,beta6,
     |  rkm7a,rkm2b,rkm2a,a1d,rkm8,rkm7b,rkm1,rkm21,rkm20,beta9,rkm24,
     |  rkm23,rkm22
      use qrj_module,only: qtotal ! qtotal(nlevp1,lon0:lon1,lat0:lat1)
      use params_module,only: spval
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: 
     |  tn,    ! neutral temperature (deg K)
     |  o2,    ! molecular oxygen (mmr)
     |  o1,    ! atomic oxygen (mmr)
     |  o1d,   ! o(1d) (mmr)
     |  o3,    ! ozone (mmr)
     |  ne,    ! electron density (cm3)
     |  no,    ! nitric oxide
     |  n4s,   ! N(4S)
     |  n2d,   ! N(2D) (updated from comp_n2d)
     |  barm,  ! mean mass
     |  xnmbar,! p0*e(-z)*barm/kT
     |  cp     ! specific heat
!
! Local:
      integer :: k,i,i0,i1
      real,parameter :: e1d=1.96
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  n2,          ! n2 (1-o2-o)/mass
     |  qnox,        ! heating from NOX chemistry
     |  qox,         ! heating from OX chemistry
     |  qo1d,        ! heating from O1D quenching and reactions
     |  qchem,       ! total heating this routine (midpoints)
     |  qchemi       ! total heating this routine (interfaces)
!
! Heating diagnostics (deg/day):
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  chemheat,    ! total heating this routine (deg/day)
     |  chmh_nox,    ! nox heating (deg/day)
     |  chmh_ox,     ! ox  heating (deg/day)
     |  chmh_o1d     ! o1d heating (deg/day)
!
      i0=lon0 ; i1=lon1
      do i=lon0,lon1
        do k=lev0,lev1-1
          n2 = (1.-o2(k,i)-o1(k,i))*rmi_n2
!
! NOX heating:
!
          qnox(k,i) = evergs*avo*xnmbar(k,i)*(n4s(k,i)*rmi_n4s*
     |     (beta1(k,i,lat)*o2(k,i)*rmi_o2*1.4+
     |      beta3(k,i,lat)*no(k,i)*rmi_no*2.68)+
     |     n2d(k,i)*rmi_n2d*(beta2*o2(k,i)*rmi_o2*1.84+
     |     beta4          *o1 (k,i)*rmi_o1*2.38+
     |     beta5(k,i,lat) *.5*(ne(k,i)+ne(k+1,i))*2.38/xnmbar(k,i)+
     |     beta6          *no (k,i)*rmi_no*5.63)+
     |     beta9 (k,i,lat)*no (k,i)*rmi_no *o3 (k,i)*rmi_o3*2.08)
!
          qchem(k,i) = qnox(k,i)
! 
! Heating from OX chemistry:
!
          qox(k,i) = evergs*avo*xnmbar(k,i)**2*
     |      (rkm20(k,i,lat)/barm(k,i)         *(o1(k,i)*rmi_o1)**2*5.12+
     |       rkm21(k,i,lat)* o2(k,i)*rmi_o2   *(o1(k,i)*rmi_o1)**2*1.10+
     |       rkm22(k,i,lat)*(o2(k,i)*rmi_o2)**2*o1(k,i)*rmi_o1    *1.10+
     |       rkm23(k,i,lat)* n2(k,i)*o2(k,i)*rmi_o2*o1(k,i)*rmi_o1*1.10+
     |       rkm24(k,i,lat)* o1(k,i)*rmi_o1    *o3(k,i)*rmi_o3    *4.06/
     |       xnmbar(k,i))
!
          qchem(k,i) = qchem(k,i)+qox(k,i)
!
! Heating from O(1D) quenching and reactions:
!
          qo1d(k,i) = evergs*avo*xnmbar(k,i)*
     |     ((rkm1 (k,i,lat)*n2 (k,i)*1.96+
     |      (rkm2a(k,i,lat)         *0.33+
     |       rkm2b(k,i,lat)         *1.96)*o2(k,i)*rmi_o2+
     |      (rkm7a*6.02+rkm7b*0.86)*o3 (k,i)*rmi_o3+
     |       rkm8*o1(k,i)*rmi_o1*1.96)*o1d(k,i)*rmi_o1)-
     |       evergs*avo*e1d*a1d*o1d(k,i)*rmi_o1
!
          qchem(k,i) = qchem(k,i)+qo1d(k,i)
!
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
      qnox (lev1,:) = spval
      qox  (lev1,:) = spval
      qo1d (lev1,:) = spval
      qchem(lev1,:) = spval
!
      call addfld('QNOX' ,' ',' ',qnox,'lev' ,lev0,lev1,'lon',i0,i1,lat)
      call addfld('QOX'  ,' ',' ',qox ,'lev' ,lev0,lev1,'lon',i0,i1,lat)
      call addfld('QO1D' ,' ',' ',qo1d,'lev' ,lev0,lev1,'lon',i0,i1,lat)
      call addfld('QCHEM',' ',' ',qchem,'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Add to total model heating:
!
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          qchemi(k,i) = qchem(k-1,i)*qchem(k,i)
          if (qchemi(k,i) < 0.) qchemi(k,i) = 1.e-20
          qtotal(k,i,lat) = qtotal(k,i,lat)+sqrt(qchemi(k,i))
        enddo ! k=lev0+1,lev1-1
!
! Lower and upper boundaries:
        qtotal(lev0,i,lat) = qtotal(lev0,i,lat)+       ! bottom
     |    1.5*qchem(lev0,i)-0.5*qchem(lev0+1,i)
        qtotal(lev1,i,lat) = qtotal(lev1,i,lat)+       ! top
     |    1.5*qchem(lev1-1,i)-0.5*qchem(lev1-2,i)
      enddo ! i=lon0,lon1
!
! Chemical heating (deg/day) for diagnostic:
      do i=lon0,lon1
        do k=lev0,lev1-1
          chemheat(k,i) = qchem(k,i)*86400./cp(k,i)
          chmh_nox(k,i) = qnox (k,i)*86400./cp(k,i)
          chmh_ox(k,i)  = qox  (k,i)*86400./cp(k,i)
          chmh_o1d(k,i) = qo1d (k,i)*86400./cp(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

      call addfld('CHEMHEAT',' ','DEG/DAY',chemheat(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('CHMH_NOX',' ','DEG/DAY',chmh_nox(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('CHMH_OX' ,' ','DEG/DAY',chmh_ox(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('CHMH_O1D',' ','DEG/DAY',chmh_o1d(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

      end subroutine qjnno
