!
      module hox_module 
!
! Advance hox by one time step. Hox family = h + oh + ho2
!
      use params_module,only: nlevp1,nlonp4,spval
      use compdat_module,only: whox ! (nlevp1)
      use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  oh_h,      ! oh/h  (ratio1)
     |  ho2_h,     ! ho2/h (ratio2)
     |  h_hox,     ! h/hox (ratio3)
     |  rmtru_hox  ! "true" mass of hox
      use addfld_module,only: addfld
      implicit none
!
      real,parameter ::
     |  phihox(3) = (/0.146, 0.243, 0.162/),
     |  alfa = -0.38,
     |  xyhox = 1.e-8
!
      real,dimension(:,:,:),allocatable :: ! (lev0:lev1,lon0:lon1,lat0:lat1)
     |  hox_psistar
      real,allocatable,dimension(:,:)   :: hox_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: hox_lbc ! lower boundary (i,3,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine alloc_hox(lon0,lon1,lat0,lat1)
!
! Allocate subdomains (without ghost cells) to module data for hox.
! This is called once per run from sub allocdata (allocdata.F).
!
      integer,intent(in) :: lon0,lon1,lat0,lat1
      integer :: istat
!
! Allocate boundary conditions:
      allocate(hox_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_hox: error allocating',
     |  ' hox_ubc: stat=',i3)") istat
      allocate(hox_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_hox: error allocating',
     |  ' hox_lbc: stat=',i3)") istat
      hox_ubc = 0. ; hox_lbc = 0.
!
! Allocate psistar:
      allocate(hox_psistar(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_hox: error allocating',
     |  ' hox_psistar: stat=',i3)") istat
      hox_psistar = 0.

      write(6,"('Completed allocations for hox.')")

      end subroutine alloc_hox
!-----------------------------------------------------------------------
      subroutine comp_hox(tn,o2,ox,op,h2o,oh,ho2,h2,h,hox,o1,o3,
     |  no,co,n4s,ch4,z,phoxic,o1d,h2o2,hplus,xnmbar,barm,hox_nm,
     |  hoxnm_out,lev0,lev1,lon0,lon1,lat)
!
! Advance HOX by one time step, and partition into ho2, oh, and h.
!
      use cons_module,only: p0,expz,expzmid,boltz,re,rmassinv_hox,
     |  rmassinv_h2o,rmassinv_h2o2,rmassinv_o1,rmassinv_h2,
     |  rmassinv_ch4,rmass_hox
! photodissociation rates from qrj are at interfaces: 
! These are interpolated to midpoints in local vars below.
      use qrj_module,only: pdh2ot,pdh2ol,pdh2o2,pdch4t
      use chemrates_module,only: rk11,rk12,rk14,rkm3,rkm4,rkm27,rkm28,
     |  rkm30,rkm31,rkm35,rkm38,rkm40,rkm41,rkm44,rkm45
      implicit none
!
! Input Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,o2,ox,op,h2o,oh,ho2,h2,h,hox,o1,o3,no,co,n4s,ch4,z,phoxic,
     |  o1d,h2o2,hplus,xnmbar,barm
      real,dimension(lev0:lev1,lon0-2:lon1+2,5),intent(in) ::
     |  hox_nm
!
! Output Args: 
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  hoxnm_out
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1
      real,dimension(lon0:lon1) :: 
     |  xnmbartop,  ! n*mbar at upper boundary (xnmbar)     (t8)
     |  hoxi,       ! hox at upper boundary interface level (t9)
     |  zkt,        ! z/Kt at top                           (t5)
     |  worki1,worki2
!
! 5/25/06 btf: changed these var names from pdxxxi to pdxxxm to
! indicate that they are on midpoints (not interfaces).
! These are interpolated from the variables use-associated
! from qrj above (they were originally calculated at interfaces by qrj)
!
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  pdh2olm,    ! pdh2ol at midpoints     (s4)
     |  pdh2otm,    ! pdh2ot at midpoints     (s5)
     |  pdh2olt,    ! pdh2otm-0.12*pdh2olm    (s6)
     |  pdh2o2m,    ! pdh2o2 at midpoints     (s14)
     |  pdch4tm,    ! pdch4t at midpoints     (s15)
     |  barmm,      ! mbar at midpoints       (s11)
     |  prod1,      ! component of hox production      (s1)
     |  prod_hox,   ! number density production of hox (s2)
     |  prod2       ! component of hox production      (s3)
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! Upper boundary condition is an upward flux of hox:
! (these are apparently not used in tgcm24)
!
      do i=lon0,lon1
        xnmbartop(i) = p0*expz(lev1-1)*expzmid*barm(lev1,i)/     ! t8
     |    (boltz*tn(lev1-1,i))
        hoxi(i) = 0.5*(hox(lev1-1,i)+hox(lev1,i))                ! t9
        zkt(i) = (re+z(lev1,i))/(9.555e4*tn(lev1-1,i))

!       write(6,"('comp_hox: i=',i3,' lat=',i3,' z(lev1,i)=',
!    |    e12.4,' tn(lev1-1,i)=',e12.4)") i,lat,z(lev1,i),tn(lev1-1,i)

        worki1(i) = 2.6183e3*hoxi(i)/rmassinv_hox*xnmbartop(i)*  ! t6
     |    sqrt(tn(lev1-1,i))*(1.+zkt(i))*exp(-zkt(i))
        worki2(i) = 2.03e-7*tn(lev1-1,i)*hoxi(i)*rmassinv_hox*   ! t7
     |    xnmbartop(i)*worki1(i)
!
! Reduce charge exchange flux to prevent H holes:
        worki2(i) = worki2(i)*0.01
        hox_ubc(i,lat) = worki1(i)+worki2(i)                     ! t4
      enddo ! i=lon0,lon1

!     write(6,"('comp_hox: lat=',i3,' xnmbartop=t8=',/,(6e12.4))") 
!    |  lat,xnmbartop
!     write(6,"('comp_hox: lat=',i3,' hoxi=t9=',/,(6e12.4))") 
!    |  lat,hoxi
!     write(6,"('comp_hox: lat=',i3,' zkt=t5=',/,(6e12.4))") 
!    |  lat,zkt
!     write(6,"('comp_hox: lat=',i3,' worki1=t6=',/,(6e12.4))") 
!    |  lat,worki1
!     write(6,"('comp_hox: lat=',i3,' worki2=t7=',/,(6e12.4))") 
!    |  lat,worki2
!     write(6,"('comp_hox: lat=',i3,' hoxubc=t4=',/,(6e12.4))") 
!    |  lat,hox_ubc(:,lat)

! Convert photodiss rates (qrj) from interfaces to midpoints:
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdh2olm(k,i) = 0.5*(pdh2ol(k,i,lat)+pdh2ol(k+1,i,lat))  ! s4
          pdh2otm(k,i) = 0.5*(pdh2ot(k,i,lat)+pdh2ot(k+1,i,lat))  ! s5
          pdh2olt(k,i) = pdh2otm(k,i)-0.12*pdh2olm(k,i)           ! s6
          pdh2o2m(k,i) = 0.5*(pdh2o2(k,i,lat)+pdh2o2(k+1,i,lat))  ! s14
          pdch4tm(k,i) = 0.5*(pdch4t(k,i,lat)+pdch4t(k+1,i,lat))  ! s15
        enddo ! k=lev0,lev1-1
! Pad top level for post-proc:
        pdh2olm(lev1,i) = spval
        pdh2otm(lev1,i) = spval
        pdh2olt(lev1,i) = spval
        pdh2o2m(lev1,i) = spval
        pdch4tm(lev1,i) = spval
      enddo ! i=lon0,lon1

!     call addfld('PDH2OLM' ,' ',' ',pdh2olm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDH2OTM' ,' ',' ',pdh2otm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDH2OLT' ,' ',' ',pdh2olt,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDH2O2M' ,' ',' ',pdh2o2m,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDCH4TM' ,' ',' ',pdch4tm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

!
! Sources:
! xnmbar == s3 == s7
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          barmm(k,i) = 0.5*(barm(k,i)+barm(k+1,i))                     ! s11
          prod1(k,i) = xnmbar(k,i)**2*((2.*rkm3*h2o(k,i)*rmassinv_h2o+ ! s1
     |      2.*rkm4*h2(k,i)*rmassinv_h2)*o1d(k,i)*rmassinv_o1+
     |      2.*rkm28(k,i,lat)*o1(k,i)*rmassinv_o1*h2(k,i)*rmassinv_h2+
     |      2.*rkm27(k,i,lat)*o1(k,i)*rmassinv_o1*h2o2(k,i)*
     |        rmassinv_h2o2)+
     |      2.*pdh2olt(k,i)*h2o (k,i)*rmassinv_h2o *xnmbar(k,i)+
     |      2.*pdh2o2m(k,i)*h2o2(k,i)*rmassinv_h2o2*xnmbar(k,i)+
     |          phoxic(k,i)*h2o (k,i)*rmassinv_h2o *xnmbar(k,i)+
     |      rk11(k,i,lat)*hplus(k,i)*o1(k,i)*rmassinv_o1*xnmbar(k,i)+
     |      pdch4tm(k,i)*ch4(k,i)*rmassinv_ch4*xnmbar(k,i)+
     |      rk14*op(k,i)*h2(k,i)*rmassinv_h2*xnmbar(k,i)+
     |      (2.*rkm44(k,i,lat)*o1(k,i)*rmassinv_o1+2.*rkm45*
     |       o1d(k,i)*rmassinv_o1)*xnmbar(k,i)*ch4(k,i)*rmassinv_ch4*
     |      xnmbar(k,i)
!
! Production of hox:
          prod_hox(k,i) = (2.*rkm38*ho2_h(k,i,lat)+2.*rkm40(k,i,lat)*  ! s2
     |      ho2_h(k,i,lat)+2.*rkm30(k,i,lat)*oh_h(k,i,lat)**2+
     |      2.*rkm31(k,i,lat)*oh_h(k,i,lat)*ho2_h(k,i,lat)+
     |      2.*rkm41(k,i,lat)*xnmbar(k,i)/barmm(k,i)+
     |      2.*rkm35(k,i,lat)*ho2_h(k,i,lat)**2)*h_hox(k,i,lat)**2
!
          prod2(k,i) = rk12*op(k,i)*h_hox(k,i,lat)                     ! s3
!
! Enforce positive:
          if (prod1   (k,i) < 1.e-20) prod1   (k,i) = 1.e-20
          if (prod_hox(k,i) < 1.e-20) prod_hox(k,i) = 1.e-20
          if (prod2   (k,i) < 1.e-20) prod2   (k,i) = 1.e-20
        enddo ! k=lev0,lev1-1
!
! Set up lower boundary such that:
!   hox_lbc(i,1,lat)*DPSX/DZ + hox_lbc(i,2,lat)*PSX + hox_lbc(i,3,lat) = 0.
! real,allocatable,dimension(:,:)   :: hox_ubc ! upper boundary (i,j)
! real,allocatable,dimension(:,:,:) :: hox_lbc ! lower boundary (i,3,j)
!
        hox_lbc(i,1,lat) = 0.                                          ! t1
        hox_lbc(i,2,lat) = 1.                                          ! t2
        worki1(i) = (-prod2(lev0,i)+sqrt(prod2(lev0,i)**2+4.*          ! t6
     |    prod_hox(lev0,i)*prod1(lev0,i)))/(2.*prod_hox(lev0,i))
        hox_lbc(i,3,lat) = -worki1(i)*rmtru_hox(lev0,i,lat)/           ! t3
     |    xnmbar(lev0,i)
      enddo ! i=lon0,lon1
!
! Only 1.e-30 diffs with tgcm24 at bottom boundary of o1d -- all others are zero diffs.
!     call addfld('XNMBAR',' ',' ',xnmbar(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2O_'  ,' ',' ',h2o   (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2_'   ,' ',' ',h2    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O1D_'  ,' ',' ',o1d   (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O1_'   ,' ',' ',o1    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2O2_' ,' ',' ',h2o2  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PHOXIC',' ',' ',phoxic(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
! H+ to lev1-1 only:
!     call addfld('HPLUS' ,' ',' ',hplus (lev0:lev1-1,i0:i1),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('CH4_'  ,' ',' ',ch4   (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
! O+ to lev1-1 only:
!     call addfld('OP_'   ,' ',' ',op    (lev0:lev1-1,i0:i1),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

!     call addfld('RKM27' ,' ',' ',rkm27 (lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RKM28' ,' ',' ',rkm28 (lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RKM44' ,' ',' ',rkm44 (lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RK11'  ,' ',' ',rk11  (lev0:lev1-1,:,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

      prod1(lev1,:) = spval
      prod2(lev1,:) = spval
      prod_hox(lev1,:) = spval
!     call addfld('PROD1'   ,' ',' ',prod1   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PROD_HOX',' ',' ',prod_hox,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PROD2'   ,' ',' ',prod2   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

!
! Do time smoothing, update hox at time n-1 (hoxnm_out), and return 
!   hox_psistar (hox at time n-1 for current step) for use in minor2:
!
      call hoxchem(barm,tn,rmtru_hox(lev0,lon0,lat),prod1,prod2,
     |  prod_hox,hox_nm,hoxnm_out,hox_psistar(lev0,lon0,lat),
     |  lev0,lev1,lon0,lon1,lat)

!     call addfld('HOXNMOUT',' ',' ',hoxnm_out(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_hox
!-----------------------------------------------------------------------
      subroutine hoxchem(barm,tn,rmtruj,prod1,prod2,prod_hox,
     |  hox_nm,hoxnm_out,psistar,lev0,lev1,lon0,lon1,lat)
!
! This subroutine deals with HOX chemistry implicitly by
!  determining psistar(HOX,n+1) from:
!
!  psistar(HOX,n+1) = (-1./(2*Dt) + sqrt(1/(2*Dt)**2+4*
!  (Q*N*MBAR/RMTRU)(n) * (psi(HOX,n-1)*dtx2inv+P*RMTRU/
!
!  From cons module:  dtx2inv = 1./(2*dt)
!
!  (N*MBAR))(n)))/(2*Q*N*MBAR/RMTRU)(n)
!
!  Here S1 = P = number density source independent of n(HOX)
!       S2 = Q = (number density loss proportional to n(HOX)**2)/n(HOX)**2
!
      use cons_module,only: 
     |  shapiro,       ! time smoothing constant
     |  p0,expz,boltz,dtx2inv
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  barm,tn
      real,dimension(lev0:lev1,lon0:lon1),intent(in) ::
     |  rmtruj,
     |  prod1,    ! s1
     |  prod_hox, ! s2
     |  prod2     ! s3
      real,dimension(lev0:lev1,lon0-2:lon1+2,5),intent(in) ::
     |  hox_nm
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  hoxnm_out
      real,dimension(lev0:lev1,lon0:lon1),intent(out) :: psistar
!
! Local:
      integer :: k,i,jm2,jm1,j,jp1,jp2,lonbeg,lonend,i0,i1,nk,nkm1
      real,dimension(lev0:lev1,lon0-2:lon1+2) :: 
     |  hoxnm1,  ! stage 1 psi(hox,n-1) (s14)
     |  hoxnm2   ! stage 2 psi(hox,n-1) (s15)
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  xnmbar,
     |  psistar1,psistar2
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
      jm2=1 ; jm1=2 ; j=3 ; jp1=4 ; jp2=5
!
! Shapiro smoother, stage 1:
      do i=lon0-2,lon1+2
!     do i=lon0,lon1
        do k=lev0,lev1
          hoxnm1(k,i) =                                       ! s14
     |      hox_nm(k,i,j)   -shapiro*(hox_nm(k,i,jp2)+
     |      hox_nm(k,i,jm2) -4.     *(hox_nm(k,i,jp1)+
     |      hox_nm(k,i,jm1))+6.*hox_nm(k,i,j))
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
! Shapiro smoother, stage 2:
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = lon1-2 ! nlonp4==nlonp4
      hoxnm2 = 0. ! init (added 5/18/04 btf)
      do i=lonbeg,lonend
        do k=lev0,lev1
          hoxnm2(k,i) =                                       ! s15
     |      hoxnm1(k,i)-shapiro*(hoxnm1(k,i+2)+hoxnm1(k,i-2)-
     |      4.*(hoxnm1(k,i+1)+hoxnm1(k,i-1))+6.*hoxnm1(k,i))
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfld('HOXNM1',' ',' ',hoxnm1(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXNM2',' ',' ',hoxnm2(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      do i=lon0,lon1
        do k=lev0,lev1-1
          xnmbar(k,i) = p0*expz(k)*0.5*(barm(k,i)+barm(k+1,i))/ ! s14
     |      (boltz*tn(k,i))
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          psistar1(k,i) = (dtx2inv+prod2(k,i))**2+4.*prod_hox(k,i)*    ! s12
     |      xnmbar(k,i)/rmtruj(k,i)*(hoxnm2(k,i)*dtx2inv+prod1(k,i)*
     |      rmtruj(k,i)/xnmbar(k,i))
          if (psistar1(k,i) < 0.) psistar1(k,i) = 1.e-20
          psistar2(k,i) = 2.*(hoxnm2(k,i)*dtx2inv+prod1(k,i)*          ! s13
     |      rmtruj(k,i)/xnmbar(k,i))/(dtx2inv+prod2(k,i)+
     |      sqrt(psistar1(k,i)))
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1
!
! Extrapolate to top level:
      do i=lon0,lon1 
        psistar2(lev1,i) = 2.*psistar2(lev1-1,i)-psistar2(lev1-2,i)
      enddo ! i=lon0,lon1 

!     call addfld('PSISTAR1',' ',' ',psistar1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('PSISTAR2',' ',' ',psistar2,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

!
! Complete previous and current hox at time n-1:
      do i=lon0,lon1
        do k=lev0,lev1
          hoxnm_out(k,i) = hox_nm(k,i,j) ! hoxnm_out is new hoxnm at itc
          psistar(k,i) = psistar2(k,i)
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfld('HOXNMITP',' ',' ',psistar           ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXNMITC',' ',' ',hoxnm_out(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine hoxchem
!-----------------------------------------------------------------------
      subroutine minor_hox(tn,o2,ox,w,difkk,hox,hox_nm,hdhox,hox_out,
     |  hoxnm_out,oh,ho2,h,lev0,lev1,lon0,lon1,lat0,lat1)
      use cons_module,only: rmass_hox,rmass_h,rmass_oh,rmass_ho2
!
! Call minor2 to complete hox. This is called from dynamics, once
! per time step after comp_hox (which is called from a latitude loop).
! Inputs are at 3d subdomains. 
!
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  o2,      ! molecular oxygen (mmr)
     |  ox,      ! oxygen family (mmr)
     |  w,       ! omega (vertical velocity)
     |  difkk,   ! eddy viscosity (from mgw)
     |  hox,     ! hox (mmr)
     |  hdhox    ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  hox_out,    ! hox output
     |  oh,         ! oh output
     |  ho2,        ! ho2 output
     |  h           ! h output
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(inout) ::
     |  hoxnm_out,  ! new hox at time n-1 (was set by hoxchem, so is intent(inout))
     |  hox_nm      ! hox at time n-1
!
! Local:
      real,parameter :: small = 1.e-12
      integer :: i,k,lat,ibndbot,ibndtop,nk,nkm1,i0,i1
      real :: hmin

      ibndbot = 0 ! ilbc in minor2
      ibndtop = 1 ! iubc in minor2
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! This is the only species that calls minor2:
      call minor2(tn,o2,ox,w,difkk,hox,hox_psistar,hdhox,hox_lbc,
     |  hox_ubc,whox,rmass_hox,phihox,alfa,xyhox,ibndbot,ibndtop,
     |  hox_out,hoxnm_out,lev0,lev1,lon0,lon1,lat0,lat1)
!
! Partition hox into ho2, oh, and h:
!     use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
!    |  oh_h,      ! oh/h  (ratio1)
!    |  ho2_h,     ! ho2/h (ratio2)
!    |  h_hox,     ! h/hox (ratio3)
!
      do lat=lat0,lat1
        do i=lon0,lon1
          do k=lev0,lev1-1

            h(k,i,lat) = h_hox(k,i,lat)*
     |        hox_out(k,i,lat)*rmass_h/rmtru_hox(k,i,lat) 

            oh(k,i,lat) = oh_h(k,i,lat)*h_hox(k,i,lat)*
     |        hox_out(k,i,lat)*rmass_oh/rmtru_hox(k,i,lat)

            ho2(k,i,lat) = ho2_h(k,i,lat)*h_hox(k,i,lat)*
     |        hox_out(k,i,lat)*rmass_ho2/rmtru_hox(k,i,lat)

            hmin = small*whox(k)
            if (h  (k,i,lat) < hmin) h  (k,i,lat) = hmin
            if (oh (k,i,lat) < hmin) oh (k,i,lat) = hmin
            if (ho2(k,i,lat) < hmin) ho2(k,i,lat) = hmin

          enddo ! k=lev0,lev1-1
!
! 10/10/06 btf: assign nlevp1 <- nlev
          h(lev1,i,lat)   = h  (lev1-1,i,lat)
          oh(lev1,i,lat)  = oh (lev1-1,i,lat)
          ho2(lev1,i,lat) = ho2(lev1-1,i,lat)
        enddo ! i=lon0,lon1

!     call addfld('HOX_OUT' ,' ',' ',hox_out(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOXNMOUT' ,' ',' ',hoxnm_out(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H_PART'  ,' ',' ',h  (:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OH_PART' ,' ',' ',oh (:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HO2_PART',' ',' ',ho2(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      enddo ! lat=lat0,lat1

      end subroutine minor_hox
!-----------------------------------------------------------------------
      end module hox_module
