!
      module co_module
      use params_module,only: nlevp1,nlonp4,nlat
      implicit none
!
!     real,dimension(nlonp4,nlat)   :: co_ubc   ! upper boundary
!     real,dimension(nlonp4,3,nlat) :: co_lbc   ! lower boundary
!     real,dimension(nlevp1,nlonp4,nlat) ::
!    |  co_prod,  ! production of co
!    |  co_loss   ! loss of co
!
! Boundary conditions and production and loss terms are allocated
! subdomains by sub alloc_co (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: co_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: co_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  co_prod,     ! production of co (k,i,j)
     |  co_loss      ! loss of co       (k,i,j)
      real :: phi_co(3) = (/0.833, 1.427, 0.852/)
      real,parameter :: alfa_co = 0. ! thermal diffusion coefficient
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_co(tn,o2,o1,barm,xnmbar,n4s,n2d,ne,o2p,oh,
     |  lev0,lev1,lon0,lon1,lat)
!
! Advance CO2 by one time step:
!
      use cons_module,only: rmass_co,p0,boltz,expz,
     |  expzmid_inv,rmassinv_o2,rmassinv_oh,rmassinv_n4s,
     |  rmassinv_o1,rmassinv_n2,rmassinv_n2d
      use chemrates_module,only: beta1,beta2,beta3,beta6,beta8,beta9,
     |  beta9n,beta17,rk5,gam14,rkm42
      use fields_module,only: tlbc ! lower boundary interface level of TN
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,         ! neutral temperature (deg K)
     |  ne,o2p,oh,
     |  n4s,n2d,
     |  o2, o1,     ! o2, o mass mixing ratios
     |  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
!
! Diffusive equilibrium at upper boundary:
      do i=lon0,lon1
	co_ubc(i,lat) = 0.
!
! no_lbc will come from comp_no.
!       no_lbc(i,1,lat) = 0.
!       no_lbc(i,2,lat) = 1.
!       no_lbc(i,3,lat) = -nob(lat)*rmass_no*boltz*tlbc(i,lat)/
!    |    (p0*expzmid_inv*expz(lev0)*barm(lev0,i))
      enddo ! i=lon0,lon1
!
! Sources:
      do i=lon0,lon1
        do k=lev0,lev1-1
          xn2(k,i) = 1.-o2(k,i)-o1(k,i)
!
	  co_prod(k,i,lat) = xnmbar(k,i)**2*o2(k,i)*rmassinv_o2*
     |      (beta1(k,i,lat)*n4s(k,i)*rmassinv_n4s+beta2*n2d(k,i)*
     |      rmassinv_n2d)
!
	  co_loss(k,i,lat) =
     |      -xnmbar(k,i)*(rkm42(k,i,lat)*oh(k,i)*rmassinv_oh+
     |      gam14(k,i,lat)*o1(k,1)*rmassinv_o1*xnmbar(k,i)/
     |      (.5*(barm(k,i)+barm(k+1,i))))

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

!     call addfsech('XNMBAR' ,' ',' ',xnmbar     ,i0,i1,nk,nkm1,lat)
!     call addfsech('CO_PROD',' ',' ',co_prod(:,i0:i1,lat),
!    |  i0,i1,nk,nkm1,lat)
!     call addfsech('CO_LOSS',' ',' ',co_loss(:,i0:i1,lat),
!    |  i0,i1,nk,nkm1,lat)
!
      end subroutine comp_co
!-----------------------------------------------------------------------
      subroutine minor_co(tn,o2,o1,co,co_nm1,co_out,co_nm1_out,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
      use cons_module,only: rmass_co
!
! Input args at full task subdomains:
      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)
     |  o1,      ! atomic oxygen (mmr)
     |  co,      ! CO  (mmr)
     |  co_nm1   ! CO  at time n-1
!
! Output args also at full task subdomains:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  co_out,    ! CO2 output
     |  co_nm1_out ! CO2 output at time n-1
!
! Local:
      integer :: lat
      integer ::nk,nkm1 ! for addfsech
!
      nk = lev1-lev0+1
      nkm1 = nk-1
!
! Minor returns co_out and co_nm1_out. Module data co_prod,
! co_loss, etc, were defined by comp_co.
!
!     subroutine minor(tn,o2,o1,fcomp,fcomp_tm1,fcomp_out,fcomp_tm1_out,
!    |  sloss,sprod,flbc,fubc,rmx,phix,alfax,lev0,lev1,lon0,lon1,lat0,
!    |  lat1,idebug)

      call minor(tn,o2,o1,co,co_nm1,co_out,co_nm1_out,
     |  co_loss,co_prod,co_lbc,co_ubc,rmass_co,phi_co,alfa_co,
     |  lev0,lev1,lon0,lon1,lat0,lat1,0)

!     do lat=lat0,lat1
!       call addfsech('CO2_OUT' ,' ',' ',co_out(:,lon0:lon1,lat),
!    |    lon0,lon1,nk,nkm1,lat)
!       call addfsech('CO2_TM1' ,' ',' ',co_nm1_out(:,lon0:lon1,lat),
!    |    lon0,lon1,nk,nkm1,lat)
!     enddo ! lat=lat0,lat1

      end subroutine minor_co
!-----------------------------------------------------------------------
      subroutine alloc_co(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(co_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_no: error allocating',
     |  ' co_ubc: stat=',i3)") istat
      allocate(co_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_no: error allocating',
     |  ' co_lbc: stat=',i3)") istat
!
! Allocate (k,i,j) subdomains to production and loss:
      allocate(co_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_no: error allocating',
     |  ' co_prod: stat=',i3)") istat
      allocate(co_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_no: error allocating',
     |  ' co_loss: stat=',i3)") istat

!     write(6,"('alloc_co: allocated module data')")
!
      end subroutine alloc_co
!-----------------------------------------------------------------------
      end module co_module