!
      module qjchem_module
!
! Chemical heating module (formerly qjnno.F).
!
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine qjchem(tn,o2,ox,no,no2,n4s,n2d,ne,o3,o1,o1d,co,co2,h2o,
     |  h2,h,oh,ho2,h2o2,ch4,cl,clo,barm,xnmbar,cp,chemheat,
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculate additions to neutral gas heating and O2 dissociation
! by NOX, HOX, OX, O(1D), CH4, and CLX chemistry. These are totalled
! in QCHEM, and added to model wide total heating QTOTAL.
!
      use qrj_module,only:   ! (nlevp1,lon0:lon1,lat0:lat1)
     |  qtotal      ! total heating                F(NQ)
      use cons_module,only:
     |  evergs,     ! 1 eV = 1.602e-12 ergs (1.602e-12)
     |  avo,        ! avogadro's number (6.023e23) 
     |  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_ho2=>rmassinv_ho2,rmi_no2=>rmassinv_no2,rmi_oh=>rmassinv_oh,
     |  rmi_o1 =>rmassinv_o1 ,rmi_h  =>rmassinv_h  ,rmi_h2=>rmassinv_h2,
     |  rmi_co =>rmassinv_co ,rmi_co2=>rmassinv_co2,
     |  rmi_ch4=>rmassinv_ch4,rmi_h2o=>rmassinv_h2o,
     |  rmi_h2o2=>rmassinv_h2o2
      use chemrates_module,only: beta1,beta2,beta3,beta4,beta5,beta6,
     |  beta8,beta9,beta10,beta11,beta12,beta13,rkm1,rkm2a,rkm2b,rkm3,
     |  rkm4,rkm5a,rkm5b,rkm6,rkm7a,rkm7b,rkm8,rkm20,rkm21,rkm22,rkm23,
     |  rkm24,rkm25,rkm26,rkm27,rkm28,rkm29,rkm30,rkm31,rkm32,rkm33,
     |  rkm34,rkm35,rkm36,rkm37,rkm38,rkm39,rkm40,rkm41,rkm42,a1d,
     |  gam1,gam2a,gam2b,gam3,gam4,gam5,gam6,gam7,gam8,gam9,gam10,
     |  gam11,gam12,gam13,gam14,gam15,del1,del2
      use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  ch3,ch3o2,ch3ooh,ch2o,cho,ch3o
      use addfld_module,only: addfld
      use params_module,only: spval
      use init_module,only: istep
      use input_module,only: step,step_compqrj
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,          ! neutral temperature
     |  o2,          ! molecular oxygen
     |  ox,          ! ox
     |  oh,          ! oh
     |  no,          ! nitric oxide
     |  no2,         ! nitric dioxide
     |  n4s,         ! n(4s)
     |  n2d,         ! n(2d)
     |  ne,          ! electron density
     |  o3, o1, o1d, ! odd oxygen
     |  co, co2,     ! carbon
     |  h2o, h2, h, ho2, h2o2, ! hydrogen, et.al.
     |  ch4,         ! methane
     |  cl,clo,      ! chlorine from solgar
     |  barm,        ! mean molecular weight   ! f(:,nj+nms)
     |  xnmbar,      ! p0e(-z)*barm/(boltz*T)  ! s1
     |  cp           ! specific heat
!
! Output args:
      real,dimension(lev0:lev1,lon0:lon1),intent(out) ::
     |  chemheat     ! total heating this routine (deg/day)
!
! Local:
      integer :: i,k,i0,i1
      real,parameter :: e1d=1.96, eo21s=0.98, eo21d=0.65, eps=0.6
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  n2,          ! n2 (1-o2-o)/mass
     |  qnox,        ! heating from NOX chemistry
     |  qhox,        ! heating from HOX chemistry
     |  qox,         ! heating from OX chemistry
     |  qo1d,        ! heating from O1D quenching and reactions
     |  qch4,        ! heating from CH4 chemistry
     |  qclx,        ! heating from CLX chemistry
     |  qchem,       ! total heating this routine (added to model qtotal) ! s2
     |  qchemi,      ! total heating this routine at interfaces
     |  barmmid      ! mean weight at midpoints
!
! Heating diagnostics (deg/day):
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  chmhnox,     ! total heating this routine (deg/day)
     |  chmhhox,     ! total heating this routine (deg/day)
     |  chmhox,      ! total heating this routine (deg/day)
     |  chmho1d,     ! total heating this routine (deg/day)
     |  chmhch4,     ! total heating this routine (deg/day)
     |  chmhclx      ! total heating this routine (deg/day)
!
      i0=lon0 ; i1=lon1
      do i=lon0,lon1 
        do k=lev0,lev1-1
          n2 = (1.-o2(k,i)-ox(k,i))*rmi_n2                  ! s3
          barmmid(k,i) = 0.5*(barm(k,i)+barm(k+1,i))
!
! NOX heating:
! (Note diffs of qnox vs tgcm24 look the same as diffs noted
!  re ox production in comp_ox.F)
!
          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*3.25)+
     |     n2d(k,i)*rmi_n2d*(beta2*o2(k,i)*rmi_o2*1.8+
     |     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)+
     |     beta8          *n4s(k,i)*rmi_n4s*oh (k,i)*rmi_oh*2.10+
     |     beta9 (k,i,lat)*no (k,i)*rmi_no *o3 (k,i)*rmi_o3*2.08+
     |     beta10(k,i,lat)*ho2(k,i)*rmi_ho2*no (k,i)*rmi_no*0.35+
     |     beta11         *o1 (k,i)*rmi_o1 *no2(k,i)*rmi_no2*1.98+
     |     beta12(k,i,lat)*o3 (k,i)*rmi_o3 *no2(k,i)*rmi_no2*1.08+
     |     beta13         *n4s(k,i)*rmi_n4s*no2(k,i)*rmi_no2*1.81)
!
          qchem(k,i) = qnox(k,i)
! 
! Heating from HOX chemistry:
!
          qhox(k,i) = evergs*avo*xnmbar(k,i)*
     |     (rkm25(k,i,lat)*o1 (k,i)*rmi_o1 *oh  (k,i)*rmi_oh  *0.72+
     |      rkm26(k,i,lat)*o1 (k,i)*rmi_o1 *ho2 (k,i)*rmi_ho2 *2.33+
     |      rkm27(k,i,lat)*o1 (k,i)*rmi_o1 *h2o2(k,i)*rmi_h2o2*3.44+
     |      rkm28(k,i,lat)*o1 (k,i)*rmi_o1 *h2  (k,i)*rmi_h2  *0.08+
     |      rkm29(k,i,lat)*oh (k,i)*rmi_oh *o3  (k,i)*rmi_o3  *1.73+
     |      rkm30(k,i,lat)*oh (k,i)*rmi_oh *oh  (k,i)*rmi_oh  *0.73+
     |      rkm31(k,i,lat)*oh (k,i)*rmi_oh *ho2 (k,i)*rmi_ho2 *3.06+
     |      rkm32(k,i,lat)*oh (k,i)*rmi_oh *h2o2(k,i)*rmi_h2o2*1.35+
     |      rkm33(k,i,lat)*oh (k,i)*rmi_oh *h2  (k,i)*rmi_h2  *0.65+
     |      rkm34(k,i,lat)*ho2(k,i)*rmi_ho2*o3  (k,i)*rmi_o3  *1.23+
     |      rkm35(k,i,lat)*1.7*(ho2(k,i)*rmi_ho2)**2+
     |      rkm36(k,i,lat)*h  (k,i)*rmi_h  *o2  (k,i)*rmi_o2*
     |        xnmbar(k,i)/barmmid(k,i)*2.11+
     |      rkm37(k,i,lat)*h  (k,i)*rmi_h  *o3  (k,i)*rmi_o3  *3.34*eps+ 
     |      rkm38         *h  (k,i)*rmi_h  *ho2 (k,i)*rmi_ho2 *2.41+
     |      rkm39(k,i,lat)*h  (k,i)*rmi_h  *ho2 (k,i)*rmi_ho2 *1.61+
     |      rkm40(k,i,lat)*h  (k,i)*rmi_h  *ho2 (k,i)*rmi_ho2 *2.34+
     |      rkm41(k,i,lat)*(h (k,i)*rmi_h)**2*
     |        xnmbar(k,i)/barmmid(k,i)*4.52+
     |      rkm42(k,i,lat)*co (k,i)*rmi_co *oh  (k,i)*rmi_oh  *1.07)
!
          qchem(k,i) = qchem(k,i)+qhox(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+
     |       rkm3          *h2o(k,i)*rmi_h2o*1.23+
     |       rkm4          *h2 (k,i)*rmi_h2 *1.88+
     |      (rkm5a*1.85+rkm5b*4.89)*ch4(k,i)*rmi_ch4+
     |       rkm6 (k,i,lat)*co2(k,i)*rmi_co2*1.96+
     |      (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)
!
! Heating from CH4 chemistry:
!
          qch4(k,i) = evergs*avo*xnmbar(k,i)*
     |      (gam1 (k,i,lat)*oh (k,i)*rmi_oh *ch4(k,i)*rmi_ch4*0.62+
     |       gam2a         *o1d(k,i)*rmi_o1 *ch4(k,i)*rmi_ch4*1.85+
     |       gam2b         *o1d(k,i)*rmi_o1 *ch4(k,i)*rmi_ch4*4.89+
! tgcm24 bug in next line used "f(i,nps1k)*rmassinv_o"
     |       gam3 (k,i,lat)*ch3(k,i,lat)*o2(k,i)*rmi_o2/
     |         barmmid(k,i)*1.22+
     |       gam15         *o1d(k,i)*rmi_o1 *co2(k,i)*rmi_co2*1.61-
     |       gam14(k,i,lat)*co (k,i)*rmi_co *ox (k,i)*rmi_o1 *
     |         xnmbar(k,i)/barmmid(k,i)*5.51-
     |       gam13(k,i,lat)*o1 (k,i)*rmi_o1 *ch4(k,i)*rmi_ch4*0.11)+
     |         evergs*avo*
     |      (gam4          *ox (k,i)*rmi_o1 *ch3(k,i,lat)*2.96+
     |       gam5 (k,i,lat)*no (k,i)*rmi_no *ch3o2(k,i,lat)*0.72+
     |       gam6 (k,i,lat)*ho2(k,i)*rmi_ho2*ch3o2(k,i,lat)*1.74+
     |       gam7 (k,i,lat)*ch3o2(k,i,lat)*ch3o2(k,i,lat)*0.24/
     |         xnmbar(k,i)+
     |       gam8 (k,i,lat)*oh (k,i)*rmi_oh*ch3ooh(k,i,lat)*1.32+ 
     |       gam9 (k,i,lat)*o2 (k,i)*rmi_o2*ch3o  (k,i,lat)*1.15+
     |       gam10         *oh (k,i)*rmi_oh*ch2o  (k,i,lat)*2.92+
     |       gam11(k,i,lat)*o1 (k,i)*rmi_o1*ch2o  (k,i,lat)*2.18-
     |       gam12(k,i,lat)*o2 (k,i)*rmi_o2*cho   (k,i,lat)*0.14)
!
          qchem(k,i) = qchem(k,i)+qch4(k,i)
!
! Heating from CLX chemistry:
!
          qclx(k,i) = evergs*avo*xnmbar(k,i)/barmmid(k,i)*
     |     (del1(k,i,lat)*cl (k,i)*o3(k,i)*rmi_o3*1.68+
     |      del2(k,i,lat)*clo(k,i)*o1(k,i)*rmi_o1*2.38)
!
          qchem(k,i) = qchem(k,i)+qclx(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1 
!
      qnox (lev1,:) = spval
      qhox (lev1,:) = spval
      qox  (lev1,:) = spval
      qo1d (lev1,:) = spval
      qch4 (lev1,:) = spval
      qclx (lev1,:) = spval
      qchem(lev1,:) = spval
!
!     call addfld('QNOX' ,' ',' ',qnox,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QHOX' ,' ',' ',qhox,'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('QCH4' ,' ',' ',qch4,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QCLX' ,' ',' ',qclx,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QCHEM',' ',' ',qchem,'lev',lev0,lev1,'lon',i0,i1,lat)
! 
! Add to total model heating:
!
! 1/2/08 btf: Increment qtotal only if qrj has been called.
!
      if (istep==1.or.mod(istep,step_compqrj/step)==0) then
        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,lev1-1
          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 
      endif

!     call addfld('QCHEMI',' ',' ',qchemi,
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QTOT_CHM',' ',' ',qtotal(:,lon0:lon1,lat),
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!
! 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)
	  chmhnox(k,i)  = qnox(k,i)*86400./cp(k,i)
	  chmhhox(k,i)  = qhox(k,i)*86400./cp(k,i)
	  chmhox(k,i)   = qox(k,i)*86400./cp(k,i)
	  chmho1d(k,i)  = qo1d(k,i)*86400./cp(k,i)
	  chmhch4(k,i)  = qch4(k,i)*86400./cp(k,i)
	  chmhclx(k,i)  = qclx(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('CHMHNOX',' ','DEG/DAY',chmhnox(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CHMHHOX',' ','DEG/DAY',chmhnox(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CHMHOX' ,' ','DEG/DAY',chmhox(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CHMHO1D',' ','DEG/DAY',chmho1d(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CHMHCH4',' ','DEG/DAY',chmhch4(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CHMHCLX',' ','DEG/DAY',chmhclx(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
      end subroutine qjchem
!-----------------------------------------------------------------------
      end module qjchem_module
