
      module o1d_module
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine comp_o1d(o2,o1,o3,n2d,ne,no,o2p,op,o1d,xnmbar,
     |  barm,lev0,lev1,lon0,lon1,lat)
!
! Update o1d.
! 9/1/06 btf: Remove o3 calculation (now using comp_major and compart,
!             where o3 is partitioned from ox)
!
      use qrj_module,only: pdo2d,pdo3d,pdo3p
      use cons_module,only: rmassinv_o2,rmassinv_o1,rmassinv_n2,
     |  rmassinv_o3,rmassinv_n2d,rmass_o1,rmassinv_no
      use chemrates_module,only: beta2,a1d,ra2,rkm2a,rkm1,rkm7a,rkm8,
     |  rkm21,rkm7b,rkm24,beta9
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
!
! Input arrays:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  o2,     ! molecular oxygen (mmr)
     |  o1,     ! atomic oxygen (mmr)
     |  o3,     ! ozone (mmr)
     |  ne,     ! electron density (cm3)
     |  no,     ! no (mmr)
     |  n2d,    ! n2d (mmr)
     |  o2p,    ! O2+ ion
     |  op,     ! O+ ion
     |  xnmbar, ! p0*e(-z)*barm/kT
     |  barm    ! mbar 
!
! Output:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  o1d     ! o(1delta) output
!
! Local:
      integer :: k,i
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pdo3dm,    ! pdo3d at midpoints
!    |  pdo3pm,    ! pdo3p at midpoints
     |  pdo2dm,    ! pdo2d at midpoints
     |  o1dprod,   ! o1d production
     |  o1dloss,   ! o1d loss
     |  o1d_cm3,   ! o1d number density
     |  o1_cm3,o2_cm3,no_cm3
      real :: psn2
!
! Photodissociation rates at midpoints:
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdo3dm(k,i) = 0.5*(pdo3d(k,i,lat)+pdo3d(k+1,i,lat))
!         pdo3pm(k,i) = 0.5*(pdo3p(k,i,lat)+pdo3p(k+1,i,lat))
          pdo2dm(k,i) = 0.5*(pdo2d(k,i,lat)+pdo2d(k+1,i,lat))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
! Scan subdomains:
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! o(1d) production:
!
	  o1dprod(k,i) = xnmbar(k,i)*
     |      (pdo2dm(k,i)*o2(k,i)*rmassinv_o2+
     |       pdo3dm(k,i)*o3(k,i)*rmassinv_o3+
     |       beta2*n2d(k,i)*rmassinv_n2d*
     |       o2(k,i)*rmassinv_o2*0.1*xnmbar(k,i))+
     |      0.85*ra2(k,i,lat)*o2p(k,i)*
     |      sqrt(ne(k,i)*ne(k+1,i))
!
! o(1d) loss:
!
	  psn2 = (1.-o1(k,i)-o2(k,i))
	  o1dloss(k,i) = xnmbar(k,i)*
     |      (rkm2a(k,i,lat) *o2 (k,i) *rmassinv_o2  +
     |       rkm1 (k,i,lat) *psn2     *rmassinv_n2  +
     |       rkm7a          *o3 (k,i) *rmassinv_o3  + 
     |       rkm8           *o1 (k,i) *rmassinv_o1) +
     |       a1d
!
! Psi o(1delta) output (assume photochemical equilibrium):
!
! Mass mixing ratio:
          o1d(k,i) = rmass_o1*o1dprod(k,i)/(o1dloss(k,i)*xnmbar(k,i))
          if (o1d(k,i) < 1.e-60) o1d(k,i) = 1.e-60
!
! Number density:
          o1d_cm3(k,i) = o1dprod(k,i)/o1dloss(k,i)
!
! o,o2,no in cm3:
          o1_cm3(k,i) = o1(k,i)*xnmbar(k,i)*rmassinv_o1
          o2_cm3(k,i) = o2(k,i)*xnmbar(k,i)*rmassinv_o2
          no_cm3(k,i) = no(k,i)*xnmbar(k,i)*rmassinv_no
!
! o3 production:
!
!         o3prod(k,i) = rkm21(k,i,lat)*o2_cm3(k,i)*o1_cm3(k,i)*
!    |      xnmbar(k,i)/barm(k,i)
!
! o3 loss:
! 8/7/06 btf: add pdo3pm
!
!         o3loss(k,i) = pdo3dm(k,i)+pdo3pm(k,i)+
!    |      (rkm7a+rkm7b)*o1d_cm3(k,i) +
!    |      rkm24(k,i,lat)*o1_cm3(k,i) +
!    |      beta9(k,i,lat)*no_cm3(k,i)
!
! Psi o3 output (assume photochemical equilibrium):
!
! Number density:
!         o3_cm3(k,i) = o3prod(k,i)/o3loss(k,i)
!
! Mass mixing ratio:
!         o3_mmr(k,i) = rmass_o3*o3_cm3(k,i)/xnmbar(k,i)
!         if (o3_mmr(k,i) < 1.e-60) o3_mmr(k,i) = 1.e-60
!         o3out(k,i) = o3_mmr(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
! 
! xnmbar:
!     call addfld('XNMBAR' ,' ',' ',xnmbar(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
! pdo3dm, pdo2dm:
!     call addfld('PDO3DM' ,' ',' ',pdo3dm(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('PDO3PM' ,' ',' ',pdo3pm(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('PDO2DM' ,' ',' ',pdo2dm(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
! o1d:
!     call addfld('O1DPROD',' ',' ',o1dprod(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O1DLOSS',' ',' ',o1dloss(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O1D'        ,' ',' ',o1d(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O1D_CM3'    ,' ',' ',o1d_cm3(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
! cm3:
!     call addfld('O2_CM3',' ',' ',o2_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O1_CM3',' ',' ',o1_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('NO_CM3',' ',' ',no_cm3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
! o3:
!     call addfld('O3PROD',' ',' ',o3prod(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O3LOSS',' ',' ',o3loss(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O3_MMR'    ,' ',' ',o3(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('O3_CM3'    ,' ',' ',o3_cm3(lev0:lev1-1,lon0:lon1),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)

      end subroutine comp_o1d
      end module o1d_module
