
      subroutine compart(o2,ox,no,o1d,o1,o3,barm,xnmbar,
     |                   o1_out,o3_out,lev0,lev1,lon0,lon1,lat)
!
! Partition OX into O and O3.
! (hydrogen is local and set to 1.0; see also comp_ox.F)
!
      use cons_module,only: rmassinv_o2,rmassinv_o1,
     |  rmassinv_h,rmassinv_no,rmass_o1,rmass_o3
      use qrj_module,only:
     |  pdo3d,  ! photodissociation of o3(d) (k,i,lat)
     |  pdo3p   ! photodissociation of o3(p) (k,i,lat)
      use chemrates_module,only: rkm21,rkm24,rkm7a,rkm7b,
     |  rkm37,beta9
      use addfld_module,only: addfld
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  o2  ,  ! updated o2 from comp_major
     |  ox  ,  ! updated ox from comp_major
     |  no  ,  ! updated no from minor_noz
     |  o1d ,  ! o(1delta) from comp_meta
     |  o1  ,  ! current o1
     |  o3  ,  ! current o3
     |  barm,  ! mbar
     |  xnmbar ! p0e^-z*mbar
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  o1_out,  ! O1 output
     |  o3_out   ! O3 output
!
! Local:
      integer :: i,k
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  o3dpsum, ! pdo3d+pdo3p
     |  o3_o,    ! o3/o ratio
     |  o_ox,    ! o/ox ratio
     |  o3_ox    ! o3/ox*(rmo1+o3/o*rmo3)
      real,dimension(lev0:lev1,lon0-2:lon1+2) ::
     |  h        ! h (h and h+ will be added to main code later, s.a. comp_ox)
      real :: temp
!
! Local hydrogen is zero for now:
      h = 0.
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          o3dpsum(k,i) = 0.5*(pdo3d(k,i,lat)+pdo3d(k+1,i,lat))+
     |                   0.5*(pdo3p(k,i,lat)+pdo3p(k+1,i,lat))
!
! O3/O ratio for partitioning of ox:
          o3_o(k,i) =
     |      xnmbar(k,i)*rkm21(k,i,lat)*o2(k,i)*rmassinv_o2*
     |      xnmbar(k,i)/(.5*(barm(k,i)+barm(k+1,i)))/(xnmbar(k,i)*
     |      (rkm24(k,i,lat)*o1 (k,i)*rmassinv_o1  +
     |      (rkm7a+rkm7b)  *o1d(k,i)*rmassinv_o1  +
     |      rkm37(k,i,lat) *h  (k,i)*rmassinv_h   +
     |      beta9(k,i,lat) *no (k,i)*rmassinv_no  +
     |      o3dpsum(k,i)))
!
! O/OX ratio: 
          o_ox(k,i) = 1./(1.+o3_o(k,i))
          if (o_ox(k,i) < 1.e-6) o_ox(k,i) = 1.e-6
          temp = (1.-o_ox(k,i))/o_ox(k,i)
          if (o3_o(k,i) < temp) o3_o(k,i) = temp
          o3_ox(k,i) = o_ox(k,i)*(rmass_o1+o3_o(k,i)*rmass_o3)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('O3DPSUM' ,' ',' ',o3dpsum(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon' ,lon0,lon1,lat)
!     call addfld('O3_O_RAT',' ',' ',o3_o(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O_OX_RAT',' ',' ',o_ox(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O3_OX'   ,' ',' ',o3_ox(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!
! Output updated o1 and o3:
      do i=lon0,lon1
        do k=lev0,lev1-1
          o1_out(k,i) = o_ox(k,i)*ox(k,i)*rmass_o1/o3_ox(k,i)
          o3_out(k,i) = o3_o(k,i)*o_ox(k,i)*ox(k,i)*rmass_o3/o3_ox(k,i)
          if (o1_out(k,i) < 1.e-10) o1_out(k,i) = 1.e-10
          if (o3_out(k,i) < 1.e-10) o3_out(k,i) = 1.e-10
        enddo ! k=lev0,lev1-1
!
! Upper boundary via log interpolation:
        o1_out(lev1,i) = o1_out(lev1-1,i)**2/o1_out(lev1-2,i)
        o3_out(lev1,i) = o3_out(lev1-1,i)**2/o3_out(lev1-2,i)
      enddo ! i=lon0,lon1

!     call addfld('O1_OUT',' ',' ',o1_out(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('O3_OUT',' ',' ',o3_out(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

      end subroutine compart
