!
      module ox_module
!
! Calculate matrix coefficients in FS for partitioning of OX into O and O3. 
! This must be called before the composition module (comp_major.F).
!
      use addfld_module,only: addfld
      use params_module,only: spval
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine comp_ox(tn,o2,ox,op,h2o,h2o2,oh,ho2,h2,h,o1,o1d,o3,no,
     |  no2,n4s,n2d,o2p,n2p,nplus,hplus,nop,ne,co,ch4,cl,clo,xiop2p,
     |  xiop2d,barm,xnmbar,lev0,lev1,lon0,lon1,lat)
!
! This routine calculates partitioning coefficients for OX in
! the array fs(lon0:lon1,nlevp1,2,0:2,lat0:lat1). The fs array
! is allocated in the chemrates module (chemrates.F).
!
      use cons_module,only: 
     |  rmassinv_o1  ,rmassinv_oh ,rmassinv_ho2,rmassinv_h  ,
     |  rmassinv_no  ,rmassinv_no2,rmassinv_n4s,rmassinv_n2d,
     |  rmassinv_h2o2,rmassinv_ch4,rmassinv_h2 ,rmassinv_co ,
     |  rmassinv_o3  ,rmassinv_o2 ,rmassinv_n2 ,
     |  rmass_o1     ,rmass_o2    ,rmass_o3
!
! Photodissociation and ionization rates from qrj, on interfaces:
      use qrj_module,only: 
     |  pdo3d,  ! photodissociation of o3(d)      (F(NJO3D))
     |  pdo3p,  ! photodissociation of o3(p)      (F(NJO3P))
     |  pdo2,   ! photodissociation of o2         (F(NRJ))
     |  qo2p,   ! o2+ ionization                  (F(NQO2P))
     |  qo2j,   !                                 (F(NQO2J)) is zero
     |  qop,    ! o+  ionization                  (F(NQOP))
     |  pdh2ol, ! lyman-alpha photodissociation of h2o  (F(NJH2OL))
     |  pdnosrb,! photodissociation of no by srb  (XJNO)
     |  pdno2   ! XJNO2
      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
     |  rmtru_ox,  ! "true" mass of ox (defined by this routine)
     |  ch3,ch2o,ch3o2,ch3o,cho
      use chemrates_module,only: beta1,beta2,beta3,beta6,beta9,beta11,
     |  beta12,beta13,rkm7a,rkm7b,rkm20,rkm21,rkm22,rkm23,rkm24,rkm25,
     |  rkm26,rkm27,rkm28,rkm29,rkm31,rkm34,rkm35,rkm36,rkm37,rkm38,
     |  rkm30,rkm40, rk1,rk3,rk4,rk5,rk6,rk7,rk8,rk9,rk10,rk11,rk12,
     |  rk16,rk23,rk24,rk27,gam3,gam4,gam6,gam7,gam9,gam11,gam12,gam13,
     |  gam14,gam15,del1,del2,ra1,ra2, 
     |  fs ! (fs(lon0:lon1,nlevp1,2,0:2,lat0:lat1) (defined by this routine)
      use init_module,only: istep ! for debug only   
      implicit none
!
! 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 from previous step
     |  op,          ! o+ ion
     |  h2o,         ! water vapor
     |  h2o2,        ! h2o2
     |  oh,          ! oh
     |  ho2,         ! ho2
     |  h2,          ! molecular hydrogen
     |  h,           ! atomic hydrogen
     |  o1,          ! atomic oxygen
     |  o1d,         ! o(1d)
     |  o3,          ! ozone
     |  no,          ! no
     |  no2,         ! no2
     |  n4s,         ! n(4s)
     |  n2d,         ! n(2d)
     |  o2p,         ! o2+ from elden
     |  n2p,         ! n2+ from elden
     |  nplus,       ! n+  from oplus
     |  hplus,       ! h+  from oplus (was xnhp)
     |  nop,         ! no+ from elden
     |  ne,          ! ne  from elden
     |  co,          ! carbon monoxide
     |  ch4,         ! methane
     |  cl,          ! from solgar
     |  clo,         ! from solgar
     |  xiop2p,      ! O+(2p) from oplus module
     |  xiop2d,      ! O+(2d) from oplus module
     |  barm,        ! mbar (interfaces)
     |  xnmbar       ! p0e(-z), etc.
!
! Local:
      integer :: i,k,i0,i1,nk,nkm1,ibndbot,ibndtop
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  o3_o,        ! o3/o ratio                 s7
     |  o_ox,        ! o/ox ratio                 s8
! 5/25/06 btf: Changed barmi to barmm to indicate midpoints
     |  barmm,       ! barm at midpoints
     |  n2,          ! n2 (mmr)                   n2=1-o2-o
! Photodissociation:
! 5/25/06 btf: Changed these from pdxxxi to pdxxxm to indicate midpoints
     |  pdh2olm,     ! pdh2ol at midpoints        s14 (NJH2OL)
     |  pdo3m,       ! pdo3d+pdo3p at midpoints   s13 (NJO3D+NJO3P)
     |  pdo2m,       ! pdo2 at midpoints          s12 (NRJ)
     |  pdno2m,      ! pdno2 at midpoints         s2  (XJNO2)
! Ionization:
     |  qo2pm,       ! qo2p at midpoints          s11 (NQO2P)
     |  qo2jm,       ! qo2j at midpoints          s10 (NQO2J)
     |  qopm,        ! qop at midpoints           s9  (NQOP)
! Production and loss rates:
     |  ox_prod1,    ! ox production              s1
     |  ox_prod2,    ! ox production              s2
     |  ox_loss1,    ! ox loss                    s3
     |  ox_loss2,    ! ox loss                    s4
     |  o2_prod1,    ! o2 production              s6
     |  o2_prod2,    ! o2 production              s9
     |  o2_prod3,    ! o2 production              s10
     |  o2_loss      ! o2 loss                    s13
!
! Diagnostics for secondary history (see addfld calls):
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  fs11,        ! fs(i,k,1,1,lat)
     |  fs12,        ! fs(i,k,1,2,lat)
     |  fs21,        ! fs(i,k,2,1,lat)
     |  fs22,        ! fs(i,k,2,2,lat)
     |  fs10,        ! fs(i,k,1,0,lat)
     |  fs20         ! fs(i,k,2,0,lat)
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! Convert several fields from full-levels (interfaces) to  half-levels 
! (midpoints). Most of these are from the qrj module.
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdh2olm(k,i)=0.5*(pdh2ol(k,i,lat)+pdh2ol(k+1,i,lat))   ! s14
          pdo3m(k,i) = 0.5*(pdo3d (k,i,lat)+pdo3d(k+1,i,lat))+   ! s13
     |                 0.5*(pdo3p (k,i,lat)+pdo3p(k+1,i,lat))
          pdo2m(k,i) = 0.5*(pdo2  (k,i,lat)+pdo2 (k+1,i,lat))    ! s12
          pdno2m(k,i)= 0.5*(pdno2 (k,i,lat)+pdno2(k+1,i,lat))    ! s2
          qo2pm(k,i) = 0.5*(qo2p  (k,i,lat)+qo2p (k+1,i,lat))    ! s11
          qo2jm(k,i) = 0.5*(qo2j  (k,i,lat)+qo2j (k+1,i,lat))    ! s10
          qopm (k,i) = 0.5*(qop   (k,i,lat)+qop  (k+1,i,lat))    ! s9
          barmm(k,i) = 0.5*(barm(k,i)+barm(k+1,i))
        enddo ! k=lev0,lev1-1
! Set top level to missing value for post-procs:
        pdh2olm(lev1,i) = spval
        pdo3m  (lev1,i) = spval
        pdo2m  (lev1,i) = spval
        pdno2m (lev1,i) = spval
        qo2pm  (lev1,i) = spval
        qo2jm  (lev1,i) = spval
        qopm   (lev1,i) = spval
      enddo ! i=lon0,lon1

! Photo-diss and ionization rates at midpoints:
!     call addfld('PDH2OLM',' ',' ',pdh2olm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDO3M'  ,' ',' ',pdo3m  ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDO2M'  ,' ',' ',pdo2m  ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('PDNO2M' ,' ',' ',pdno2m ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QO2PM'  ,' ',' ',qo2pm  ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QO2JM'  ,' ',' ',qo2jm  ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('QOPM'   ,' ',' ',qopm   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      do i=lon0,lon1
        do k=lev0,lev1-1
!     
! Ratio O3/O:
!
          o3_o(k,i) = xnmbar(k,i)*rkm21(k,i,lat)*o2(k,i)*rmassinv_o2* ! s7
     |      xnmbar(k,i)/barmm(k,i) / (xnmbar(k,i)*
     |      (rkm24(k,i,lat)*o1 (k,i)*rmassinv_o1+
     |      (rkm7a+rkm7b)  *o1d(k,i)*rmassinv_o1+
     |      rkm29(k,i,lat) *oh (k,i)*rmassinv_oh+
     |      rkm34(k,i,lat) *ho2(k,i)*rmassinv_ho2+
     |      rkm37(k,i,lat) *h  (k,i)*rmassinv_h+
     |      beta9(k,i,lat) *no (k,i)*rmassinv_no+
     |      beta12(k,i,lat)*no2(k,i)*rmassinv_no2)+
     |      pdo3m(k,i)+del1(k,i,lat)*cl(k,i)*xnmbar(k,i)/barmm(k,i))
!
! Ratio O/OX:
          o_ox(k,i) = 1./(1.+o3_o(k,i))                        ! s8
!
! OX production terms:
!
          n2(k,i) = 1.-o2(k,i)-ox(k,i)
          ox_prod1(k,i) = xnmbar(k,i)**2*                      ! s1
     |      (rkm30(k,i,lat)*(oh(k,i)*rmassinv_oh)**2+
     |       rkm40(k,i,lat)*h  (k,i)*rmassinv_h  *ho2(k,i)*rmassinv_ho2+
     |       beta3(k,i,lat)*n4s(k,i)*rmassinv_n4s*no (k,i)*rmassinv_no +
     |       beta6         *n2d(k,i)*rmassinv_n2d*no (k,i)*rmassinv_no +
     |       beta13        *n4s(k,i)*rmassinv_n4s*no2(k,i)*rmassinv_no2)
     |      +0.5*(pdnosrb(k,i,lat)+pdnosrb(k+1,i,lat))*
     |       no(k,i)*rmassinv_no*xnmbar(k,i)+
     |       pdno2m(k,i)*no2(k,i)*rmassinv_no2*xnmbar(k,i)+
     |       xnmbar(k,i)*(rk4*o2p(k,i)*n4s(k,i)*rmassinv_n4s 
     |      +rk10*op(k,i)*n2d(k,i)*rmassinv_n2d+
     |       rk12*op(k,i)*h(k,i)*rmassinv_h+
     |       (rk16*xiop2p(k,i)+rk24*xiop2d(k,i))*n2(k,i)*rmassinv_n2)+
     |       (ra1(k,i,lat)*nop(k,i)+2.*ra2(k,i,lat)*o2p(k,i))*
     |       sqrt(ne(k,i)*ne(k+1,i))
!
          ox_prod2(k,i) = xnmbar(k,i)*                         ! s2
     |     (beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+
     |      beta2*n2d(k,i)*rmassinv_n2d)+
     |      rk1(k,i,lat)*op(k,i)+rk7*nplus(k,i)+      
     |      rk23*xiop2p(k,i)+rk27*xiop2d(k,i)+
     |      2.*pdo2m(k,i)
!
! OX loss terms:
!
          ox_loss1(k,i) =                                      ! s3
     |      (2.*rkm20(k,i,lat)*xnmbar(k,i)/barmm(k,i)+
     |       2.*rkm24(k,i,lat)*o3_o(k,i))*o_ox(k,i)**2
!
          ox_loss2(k,i) = (xnmbar(k,i)*                        ! s4
     |     (rkm29(k,i,lat)*oh  (k,i)*rmassinv_oh *o3_o(k,i)+
     |      rkm34(k,i,lat)*ho2 (k,i)*rmassinv_ho2*o3_o(k,i)+
     |      rkm37(k,i,lat)*h   (k,i)*rmassinv_h  *o3_o(k,i)+       
     |      rkm25(k,i,lat)*oh  (k,i)*rmassinv_oh  +
     |      rkm26(k,i,lat)*ho2 (k,i)*rmassinv_ho2 +
     |      rkm27(k,i,lat)*h2o2(k,i)*rmassinv_h2o2+
     |      gam13(k,i,lat)*ch4 (k,i)*rmassinv_ch4 +
     |      rkm28(k,i,lat)*h2  (k,i)*rmassinv_h2  +
     |      gam14(k,i,lat)*xnmbar(k,i)/barmm(k,i)*co(k,i)*rmassinv_co+
     |      beta9 (k,i,lat)*no (k,i)*rmassinv_no *o3_o(k,i)+
     |      beta12(k,i,lat)*no2(k,i)*rmassinv_no2*o3_o(k,i)+
     |      beta11         *no2(k,i)*rmassinv_no2)+
     |      gam4        *ch3(k,i,lat)+gam11(k,i,lat)*ch2o(k,i,lat)+
     |      rk3(k,i,lat)*n2p(k,i)    +rk8           *nplus(k,i)+
     |      rk11(k,i,lat)*hplus(k,i))*o_ox(k,i)+
     |     (del1(k,i,lat)*cl   (k,i) *o3_o(k,i)*o_ox(k,i)+
     |      del2(k,i,lat)*clo  (k,i) *o_ox(k,i))*
     |      xnmbar(k,i)/barmm(k,i)
!
! O2 production:
!
          o2_prod1(k,i) =                                      ! s6
     |      (rkm20(k,i,lat)*xnmbar(k,i)/barmm(k,i)+
     |      2.*rkm24(k,i,lat)*o3_o(k,i))*o_ox(k,i)**2
!
          o2_prod2(k,i) =                                      ! s9
     |      ((rkm25(k,i,lat)*oh (k,i)*rmassinv_oh+
     |        rkm26(k,i,lat)*ho2(k,i)*rmassinv_ho2+
     |        rkm29(k,i,lat)*oh (k,i)*rmassinv_oh *o3_o(k,i)+
     |     2.*rkm34(k,i,lat)*ho2(k,i)*rmassinv_ho2*o3_o(k,i)+
     |        rkm37(k,i,lat)*h  (k,i)*rmassinv_h  *o3_o(k,i)+
     |        beta9(k,i,lat)*no (k,i)*rmassinv_no *o3_o(k,i)+
     |       beta11         *no2(k,i)*rmassinv_no2+
     |       beta12(k,i,lat)*no2(k,i)*rmassinv_no2*o3_o(k,i))*
     |         xnmbar(k,i)+pdo3m(k,i)*o3_o(k,i))*o_ox(k,i)+
     |       (del1(k,i,lat) *cl (k,i)*o3_o(k,i) *o_ox(k,i)+
     |        del2(k,i,lat) *clo(k,i)*o_ox(k,i))*xnmbar(k,i)/barmm(k,i)
!
          o2_prod3(k,i) = xnmbar(k,i)**2*                      ! s10
     |      (rkm31(k,i,lat)*oh(k,i)*rmassinv_oh*ho2(k,i)*rmassinv_ho2+
     |       rkm38         *h (k,i)*rmassinv_h *ho2(k,i)*rmassinv_ho2+
     |       rkm35(k,i,lat)*(ho2(k,i)*rmassinv_ho2)**2)+
     |       rk5*no(k,i)*rmassinv_no*o2p(k,i)*xnmbar(k,i)+
     |       gam6(k,i,lat)*ch3o2(k,i,lat)*ho2(k,i)*rmassinv_ho2*
     |         xnmbar(k,i)+gam7(k,i,lat)*ch3o2(k,i,lat)**2+
     |       gam15*co(k,i)*rmassinv_co*o1d(k,i)*rmassinv_o1*
     |         xnmbar(k,i)**2
!
! O2 loss:
!
          o2_loss(k,i) = xnmbar(k,i)*                          ! s13
     |      (rkm36(k,i,lat)*xnmbar(k,i)/barmm(k,i)*h(k,i)*rmassinv_h+
     |       beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+
     |       beta2         *n2d(k,i)*rmassinv_n2d)+
     |       rk1(k,i,lat)*op(k,i)+(rk6+rk7)*nplus(k,i)+
     |       rk9*n2p(k,i)+pdo2m(k,i)+rk23*xiop2p(k,i)+rk27*xiop2d(k,i)+
     |       gam3(k,i,lat)*ch3(k,i,lat)*xnmbar(k,i)/barmm(k,i)+
     |       gam9(k,i,lat)*ch3o(k,i,lat)+gam12(k,i,lat)*cho(k,i,lat) 
!
! "True" OX mass:
!
          rmtru_ox(k,i,lat) = o_ox(k,i)*(rmass_o1+o3_o(k,i)*rmass_o3) ! s7
!
! Matrix coefficients (fs(lon0:lon1,nlevp1,2,0:2,lat0:lat1))
!
          fs(i,k,1,1,lat) = -o2_loss(k,i)-xnmbar(k,i)**2*o_ox(k,i)*
     |      ox(k,i)/rmtru_ox(k,i,lat)*(1./3.*rkm21(k,i,lat)*o_ox(k,i)*
     |      ox(k,i)/rmtru_ox(k,i,lat)+ 2./3.*rkm22(k,i,lat)*
     |      o2(k,i)*rmassinv_o2      + 1./2.*rkm23(k,i,lat)*
     |      n2(k,i)*rmassinv_n2)
          fs(i,k,1,2,lat) = (o2_prod2(k,i)+xnmbar(k,i)*o2_prod1(k,i)*
     |      ox(k,i)/rmtru_ox(k,i,lat)-xnmbar(k,i)**2*o_ox(k,i)*
     |      o2(k,i)*rmassinv_o2*(2./3.*rkm21(k,i,lat)*o_ox(k,i)*
     |      ox(k,i)/rmtru_ox(k,i,lat)+1./3.*rkm22(k,i,lat)*o2(k,i)*
     |      rmassinv_o2+1./2.*rkm23(k,i,lat)*n2(k,i)*rmassinv_n2))*
     |      rmass_o2/rmtru_ox(k,i,lat)
          fs(i,k,2,1,lat) =  ox_prod2(k,i)*rmtru_ox(k,i,lat)*rmassinv_o2
          fs(i,k,2,2,lat) = -ox_loss2(k,i)-ox_loss1(k,i)*ox(k,i)/
     |      rmtru_ox(k,i,lat)*xnmbar(k,i)
          fs(i,k,1,0,lat) = (o2_prod3(k,i)-qo2pm(k,i))*rmass_o2/
     |      xnmbar(k,i)
          fs(i,k,2,0,lat) = (ox_prod1(k,i)-qopm(k,i))*rmtru_ox(k,i,lat)/
     |      xnmbar(k,i)
!
! Save fs to secondary history:
          fs11(k,i) = fs(i,k,1,1,lat)
          fs12(k,i) = fs(i,k,1,2,lat)
          fs21(k,i) = fs(i,k,2,1,lat)
          fs22(k,i) = fs(i,k,2,2,lat)
          fs10(k,i) = fs(i,k,1,0,lat)
          fs20(k,i) = fs(i,k,2,0,lat)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Save diagnostics up to lev1-1:
!
!     call addfld('O3_O'    ,' ',' ',o3_o(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('OX_PROD1',' ',' ',ox_prod1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('OX_PROD2',' ',' ',ox_prod2(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('OX_LOSS1',' ',' ',ox_loss1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('OX_LOSS2',' ',' ',ox_loss2(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O2_PROD1',' ',' ',o2_prod1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O2_PROD2',' ',' ',o2_prod2(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O2_PROD3',' ',' ',o2_prod3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O2_LOSS' ,' ',' ',o2_loss(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RMTRU_OX',' ',' ',
!    |  rmtru_ox(lev0:lev1-1,lon0:lon1,lat),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
      fs11(lev1,:) = spval ; fs12(lev1,:) = spval
      fs21(lev1,:) = spval ; fs22(lev1,:) = spval
      fs10(lev1,:) = spval ; fs20(lev1,:) = spval
!     call addfld('FS11' ,' ',' ',fs11,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS12' ,' ',' ',fs12,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS12' ,' ',' ',fs12,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS21' ,' ',' ',fs21,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS21' ,' ',' ',fs21,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS22' ,' ',' ',fs22,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS22' ,' ',' ',fs22,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS10' ,' ',' ',fs10,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS10' ,' ',' ',fs10,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS20' ,' ',' ',fs20,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FS20' ,' ',' ',fs20,'lev',lev0,lev1,'lon',i0,i1,lat)
 
      end subroutine comp_ox
!-----------------------------------------------------------------------
      end module ox_module