!
      subroutine hoxpion(tn,o2,ox,o1,no,h2o,co2,ne,barm,xnmbar,
     |  phoxic,solar_protons_flag,lev0,lev1,lon0,lon1,lat)
!
! HOX production from D-region ion chemistry.
! This is calculated only if solar_protons_flag (namelist read 
!   parameter solar_protons) is > 0. (default is solar_protons==0)
!
      use qrj_module,only: qo2p,qop,qn2p,qnop,qnp
      use params_module,only: zibot,dlev,spval
      use chemrates_module,only: rin40,ra2,rk5,rin9,rin10,rin11,
     |  rin12,rin13,rin17,rin18,rin19,rin20,rin21,ra1
      use cons_module,only: rmassinv_no,rmassinv_o2,rmassinv_h2o,
     |  rmassinv_o1,rmassinv_n2,rmassinv_co2
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: solar_protons_flag,lev0,lev1,
     |  lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,o2,ox,o1,no,h2o,co2,ne,barm,xnmbar
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  phoxic
!
! Local:
      real,parameter :: zptop = -7. ! top of phoxic calculation
      real,parameter ::
     |  rin14=1.E-9, rin15=1.E-9  ,rin16=1.E-9  ,rin22=1.e-9  ,
     |  rin23=1.E-9, rin24=1.E-9  ,rin25=7.E-11 ,rin41=1.5e-9 ,
     |  rin44=1.E-9, rin45=2.E-10 ,rin48=3.E-10 ,rin49=4.4e-10,
     |  alph2=3.E-6
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  work ,work1,work2,work3,work4,work5,work6,work7,work8,
     |  work9,work10,work11,work12,work13,work14,work15,work16,
     |  barmi,
     |  prod_ion  ! total production rate of ion species 
      integer :: k,i,ktop,i0,i1,nk,nkm1
      real :: zp,zp1,psn2
!
      i0 = lon0 ; i1 = lon1 ; nk = lev1-lev0+1 ; nkm1 = nk-1
!
! Do not calculate phoxic if not doing solar protons:
!     if (solar_protons_flag <= 0) then
!       phoxic(:,:) = 1.e-60
!       prod_ion(:,:) = 1.e-60
!       return
!     endif
!
! Summation over ion species production rates:
      work = 0. ! whole-array init                     ! s1
      work(:,:) = work(:,:)+qo2p(:,:,lat)
      work(:,:) = work(:,:)+qop (:,:,lat)
      work(:,:) = work(:,:)+qn2p(:,:,lat)
      work(:,:) = work(:,:)+qnop(:,:,lat)
      work(:,:) = work(:,:)+qnp (:,:,lat)
!
! Interfaces:
      do i=lon0,lon1
        do k=lev0,lev1-1
          prod_ion(k,i) = 0.5*(work(k,i)+work(k+1,i))  ! s2
          barmi(k,i) = 0.5*(barm(k,i)+barm(k+1,i))     ! s9
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Find index to top zp level:
      ktop = 0
      do k=lev0,lev1-1
        zp  = zibot+(k-1)*dlev
        zp1 = zibot+(k  )*dlev
        if (zptop >= zp .and. zptop <= zp1) ktop = k
      enddo
      if (ktop==0) then
        write(6,"(/,'>>> hoxpion: cannot find k-index to zptop=',
     |    f8.2)") zptop
        call shutdown('hoxpion')
      endif
!
! These inits are not really necessary, but help when comparing 
! with tgcm24 because these are defined only up to zptop.
      work1 = 0. ; work2 = 0. ; work3 = 0. ; work4 = 0.
      work5 = 0. ; work6 = 0. ; work7 = 0. ; work8 = 0.
      work9 = 0. ; work10 =0. ; work11= 0. ; work12= 0.
      work13= 0. ; work14 =0. ; work15= 0. ; work16= 0.
      phoxic = 0.
!
      do i=lon0,lon1
        do k=lev0,ktop
!
! Ion production of hox from O2+ ion chemical chain:
!
          work1(k,i) =rin40(k,i,lat)*o2(k,i)*rmassinv_o2*         ! s10
     |      xnmbar(k,i)**2/barmi(k,i)/(ra2(k,i,lat)*ne(k,i)+
     |      rin40(k,i,lat)*o2(k,i)*rmassinv_o2*xnmbar(k,i)**2/
     |      barmi(k,i)+rk5*no(k,i)*rmassinv_no*xnmbar(k,i))
!
          work2(k,i) = rin41/(rin48*o1(k,i)*rmassinv_o1*xnmbar(k,i)+  ! s11
     |      alph2*ne(k,i)+rin41*h2o(k,i)*rmassinv_h2o*xnmbar(k,i))
!
          work3(k,i) = 
     |       (rin44+rin45)*h2o(k,i)*rmassinv_h2o*xnmbar(k,i)/     ! s12
     |      ((rin44+rin45)*h2o(k,i)*rmassinv_h2o*xnmbar(k,i)+
     |      alph2*ne(k,i))
!
          work4(k,i) = 2.*                                        ! s13
     |      (0.904+rin48*o1(k,i)*rmassinv_o1*xnmbar(k,i)*
     |      (0.904*work1(k,i)/(alph2*ne(k,i)+
     |      rin41*h2o(k,i)*rmassinv_h2o*xnmbar(k,i)+
     |      (1.-work1(k,i))*rin48*o1(k,i)*rmassinv_o1*xnmbar(k,i))))*
     |      work1(k,i)*work2(k,i)*work3(k,i)
!
! Ion production of hox from NO+ ion chemical chain:
!
          psn2 = 1.-o2(k,i)-ox(k,i)
          work5(k,i) =                                            ! s1
     |      (rin9 (k,i,lat)*psn2    *rmassinv_n2+
     |       rin11(k,i,lat)*co2(k,i)*rmassinv_co2+
     |       rin13(k,i,lat)*h2o(k,i)*rmassinv_h2o)*
     |      xnmbar(k,i)**2/barmi(k,i)+alph2*ne(k,i)
!
          work6(k,i) =                                            ! s4
     |      rin12(k,i,lat)*xnmbar(k,i)/barmi(k,i)+
     |      rin15*h2o(k,i)*rmassinv_h2o*xnmbar(k,i)+
     |      alph2*ne(k,i)
!
          work7(k,i) =                                            ! s5
     |       rin10(k,i,lat)*xnmbar(k,i)/barmi(k,i)+
     |      (rin14*co2(k,i)*rmassinv_co2+
     |       rin16*h2o(k,i)*rmassinv_h2o)*xnmbar(k,i)+
     |       alph2*ne(k,i)
!
          work8(k,i) =                                            ! s6
     |      rin9(k,i,lat)*psn2*rmassinv_n2*xnmbar(k,i)**2/barmi(k,i)
!
          work9(k,i) =                                           ! s7
     |      1.-(alph2*ne(k,i)/(work5(k,i)-xnmbar(k,i)/barmi(k,i)*
     |      (rin10(k,i,lat)*work8(k,i)/work7(k,i)+rin12(k,i,lat)*
     |      (rin11(k,i,lat)*co2(k,i)*rmassinv_co2*xnmbar(k,i)**2/
     |      barmi(k,i)+work8(k,i)*rin14*co2(k,i)*rmassinv_co2*
     |      xnmbar(k,i))/work7(k,i)/work6(k,i))))
!
          work10(k,i) =                                          ! s8
     |      (rin17(k,i,lat)*psn2    *rmassinv_n2+
     |       rin19(k,i,lat)*co2(k,i)*rmassinv_co2+
     |       rin21(k,i,lat)*h2o(k,i)*rmassinv_h2o)*
     |      xnmbar(k,i)**2/barmi(k,i)+alph2*ne(k,i)
!
          work11(k,i) =                                          ! s14
     |      rin23*h2o(k,i)*rmassinv_h2o*xnmbar(k,i)+
     |      rin20(k,i,lat)*xnmbar(k,i)/barmi(k,i)+alph2*ne(k,i)
!
          work12(k,i) =                                          ! s15
     |      (rin22*co2(k,i)*rmassinv_co2+
     |       rin24*h2o(k,i)*rmassinv_h2o)*xnmbar(k,i)+
     |       rin18(k,i,lat)*xnmbar(k,i)/barmi(k,i)+alph2*ne(k,i)
!
          work13(k,i) =                                          ! s1
     |      rin17(k,i,lat)*psn2*rmassinv_n2*xnmbar(k,i)**2/barmi(k,i)
!
          work14(k,i) =                                          ! s4
     |      1.-(alph2*ne(k,i)/(work10(k,i)-xnmbar(k,i)/barmi(k,i)*
     |      (rin18(k,i,lat)*work13(k,i)/work12(k,i)+
     |       rin20(k,i,lat)*(rin19(k,i,lat)*co2(k,i)*rmassinv_co2*
     |       xnmbar(k,i)**2/barmi(k,i)+work13(k,i)*rin22*co2(k,i)*
     |       rmassinv_co2*xnmbar(k,i))/work12(k,i)/work11(k,i))))
!
          work15(k,i) = rin25/(rin25*h2o(k,i)*rmassinv_h2o*      ! s5
     |      xnmbar(k,i)+alph2*ne(k,i))
!
          work16(k,i) = 2.*work9(k,i)*work14(k,i)*work15(k,i)*   ! s10
     |      (.096+.904*rin49*xnmbar(k,i)*no(k,i)*rmassinv_no*
     |      (1.+ra1(k,i,lat)*rin48*o1(k,i)*rmassinv_o1*
     |      xnmbar(k,i)/(alph2*ne(k,i)+rin41*h2o(k,i)*rmassinv_h2o*
     |      xnmbar(k,i)+(1.-ra1(k,i,lat))*rin48*o1(k,i)*rmassinv_o1*
     |      xnmbar(k,i)))/(alph2*ne(k,i)+rin40(k,i,lat)*o2(k,i)*
     |      rmassinv_o2*xnmbar(k,i)**2/barmi(k,i)+rin49*no(k,i)*
     |      rmassinv_no*xnmbar(k,i)))
!
! Total hox production output:
          phoxic(k,i) = (work4(k,i)+work16(k,i))*prod_ion(k,i)
        enddo ! k=lev0,ktop
      enddo ! i=lon0,lon1
!
! These addfld calls will executed only if solar_protons is set,
! see return statement above. 
!

      call addfld('PHOXNO',' ',' ',work1,
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK2',' ',' ',work2,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK3',' ',' ',work3,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('PHOXO2',' ',' ',work4,
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK5',' ',' ',work5,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK6',' ',' ',work6,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK7',' ',' ',work7,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK8',' ',' ',work8,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORK9',' ',' ',work9,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKa',' ',' ',work10,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKb',' ',' ',work11,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKc',' ',' ',work12,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKd',' ',' ',work13,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKe',' ',' ',work14,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKf',' ',' ',work15,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXWORKg',' ',' ',work16,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('PHOXIC'  ,' ',' ',phoxic(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
      prod_ion(lev1,:) = spval
      call addfld('QTOT_AUR'  ,' ',' ',prod_ion(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine hoxpion
