!
      module o1d_module
      use params_module,only: nlevp1,nlonp4,nlat
      implicit none
!
!     real,dimension(nlonp4,nlat)   :: o1d_ubc   ! upper boundary
!     real,dimension(nlonp4,3,nlat) :: o1d_lbc   ! lower boundary
!     real,dimension(nlevp1,nlonp4,nlat) ::
!    |  o1d_prod,  ! production of o1d
!    |  o1d_loss   ! loss of o1d
!
!
      real,allocatable,dimension(:,:,:) ::
     |  o1d_prod,     ! production of o1d (k,i,j)
     |  o1d_loss      ! loss of o1d       (k,i,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_o1d(tn,o2,o1,n4s,n2d,ne,o2p,op,o1d,co,oh,
     |  barm,xnmbar,lev0,lev1,lon0,lon1,lat)
!
! Advance O1D by one time step:
!
      use cons_module,only: rmass_o1d,p0,boltz,expz,
     |  expzmid_inv,rmassinv_o2,rmassinv_co,rmassinv_oh,
     |  rmassinv_o1,rmassinv_n2
      use chemrates_module,only: beta1,beta2,beta3,beta6,beta8,beta9,
     |  beta9n,beta17,rk5,rk13,gam14,rkm42,gam15,o1dmix
      use fields_module,only: tlbc ! lower boundary interface level of TN
      use qrj_module,only:
     |  pdo1dt    ! total photodissociation of o1d
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,           ! neutral temperature (deg K)
     |  o2, o1,       ! o2, o mass mixing ratios
     |  n4s, n2d,     ! n4s, n2d mass mixing ratios
     |  ne, o2p,      ! electron density, o2+
     |  op,o1d,co,oh, ! additional inputs
     |  barm,         ! mean molecular weight
     |  xnmbar        ! p0*e(-z)*barm/kT
!
! Local:
      integer :: k,i
      integer ::i0,i1,nk,nkm1 ! for addfsech
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  xn2       ! N2 (mmr)
!
      i0 = lon0
      i1 = lon1
      nk = lev1-lev0+1
      nkm1 = nk-1
!
!
! Sources:
      do i=lon0,lon1
        do k=lev0,lev1-1
          xn2(k,i) = 1.-o2(k,i)-o1(k,i)
!



!
! Metastable for o(1d)
      do i=lon0,lon1
        do k=lev0,lev1-1
	  pdo3di(k,i,lat) = 0.5*(pdo3d(k,i,lat)+pdo3d(k+1,i,lat))
	  pdo2di(k,i,lat) = 0.5*(pdo2d(k,i,lat)+pdo2d(k+1,i,lat))
!         pdco2di(k,i,lat) = 0.5*(pdco2d(k,i,lat)+pdco2d(k+1,i,lat))
!         pdco2di(k,i,lat) = 0.
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1


      do i=lon0,lon1
        do k=lev0,lev1-1
!
! o(1d) production:
!
	  o1d_prod(k,i,lat) = xnmbar(k,i)*
     |      (pdo2di(k,i,lat)*o2 (k,i)*rmassinv_o2+
     |       pdo3di(k,i,lat)*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))
	  o1d_loss(k,i,lat) = 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  +
     |       gam15          *co2(k,i) *rmassinv_co2)+
     |       a1d
!
! Psi o(1delta) output (was xno1d in tgcm24):
!
	  o1d(k,i) = rmass_o1*o1d_prod(k,i,lat)/(o1d_loss(k,i)*
     |               xnmbar(k,i))
          if (o1d(k,i) < 1.e-60) o1d(k,i) = 1.e-60
!

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfsech('XNMBAR' ,' ',' ',xnmbar     ,i0,i1,nk,nkm1,lat)
!     call addfsech('O1D_PROD',' ',' ',o1d_prod(:,i0:i1,lat),
!    |  i0,i1,nk,nkm1,lat)
!     call addfsech('O1D_LOSS',' ',' ',o1d_loss(:,i0:i1,lat),
!    |  i0,i1,nk,nkm1,lat)
!
      end subroutine comp_o1d
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      subroutine alloc_o1d(lon0,lon1,lat0,lat1)
!
! Allocate subdomains (without ghost cells) to module data for boundary 
! conditions and production and loss terms. This is called once per run 
! from sub allocdata (allocdata.F).
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate (i,j) subdomains to boundary conditions:
      allocate(o1d_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_o1d: error allocating',
     |  ' o1d_ubc: stat=',i3)") istat
      allocate(o1d_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_o1d: error allocating',
     |  ' o1d_lbc: stat=',i3)") istat
!
! Allocate (k,i,j) subdomains to production and loss:
      allocate(o1d_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_o1d: error allocating',
     |  ' o1d_prod: stat=',i3)") istat
      allocate(o1d_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_o1d: error allocating',
     |  ' o1d_loss: stat=',i3)") istat

!     write(6,"('alloc_o1d: allocated module data')")
!
      end subroutine alloc_o1d
!-----------------------------------------------------------------------
      end module o1d_module
