!
      module co_module
!
! Advance co by one time step.
!
! Boundary conditions, production and loss for CO are defined
! by comp_co, and referenced by minor_co. Comp_co is called
! from a latitude loop in dynamics. After comp_co, dynamics calls
! minor_co, which passes this module data to sub minor. Sub
! minor contains 3d mpi calls and its own latitude loops. 
!
      use params_module,only: nlevp1
      use addfld_module,only: addfld
      implicit none
!
! 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)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_co(o1,co,co2,hox,o1d,op,ch4,cl,barm,xnmbar,
     |  lev0,lev1,lon0,lon1,lat)
      use solgar_module,only: xcolb ! (nlat)
      use cons_module,only: rmass_co,rmassinv_o1,rmassinv_co,
     |  rmassinv_co2,rmassinv_ch4
      use qrj_module,only:
     |  pdco2t    ! total photodissociation of co2  (F(NJCO2T))
      use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  oh_h,     ! oh/h  (ratio1)
     |  h_hox,    ! h/hox (ratio3)
     |  rmtru_hox ! "true" mass of hox
      use chemrates_module,only: rk13,rkm5a,rkm5b,rkm42,rkm43,rkm44,
     |  rkm45,gam14,gam15,del3
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  o1,          ! atomic oxygen
     |  co,          ! carbon monoxide
     |  co2,         ! carbon dioxide
     |  hox,         ! hox
     |  o1d,         ! o1(delta)
     |  op,          ! o+
     |  ch4,         ! methane 
     |  cl,          ! cl from comp_solgar
     |  barm,        ! mbar
     |  xnmbar       ! p0e(-z), etc.
!
! Local:
      integer :: i,k,i0,i1,nk,nkm1,ibndbot,ibndtop
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pdco2ti      ! photodissociation of co2 at midpoints
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
      do i=lon0,lon1 
        co_lbc(i,1,lat) = 0.    ! t1
        co_lbc(i,2,lat) = 1.    ! t2
        co_lbc(i,3,lat) = -xcolb(lat)*rmass_co/barm(lev0,i)
        co_ubc(i,lat)   = 0.
      enddo
!     write(6,"('comp_co: lat=',i3,' co_lbc(3)=',/,(6e12.4))") 
!    |  lat,co_lbc(:,3,lat)
!
! Photodiss at midpoints:
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdco2ti(k,i) = 0.5*(pdco2t(k,i,lat)+pdco2t(k+1,i,lat)) ! s4
!
! Number density production of co:
!
          co_prod(k,i,lat) = xnmbar(k,i)*(pdco2ti(k,i)+rk13*op(k,i)+
     |      gam15*o1d(k,i)*rmassinv_o1*xnmbar(k,i))*co2(k,i)*
     |      rmassinv_co2+(rkm43(k,i,lat)*hox(k,i)/rmtru_hox(k,i,lat)*
     |      oh_h(k,i,lat)*h_hox(k,i,lat)+rkm45*o1d(k,i)*rmassinv_o1+
     |      rkm44(k,i,lat)*o1(k,i)*rmassinv_o1+del3(k,i,lat)*cl(k,i)/
     |      (.5*(barm(k,i)+barm(k+1,i))))*ch4(k,i)*rmassinv_ch4*
     |      xnmbar(k,i)*xnmbar(k,i)
!
! Number density loss of co:
!
          co_loss(k,i,lat) = -xnmbar(k,i)*(rkm42(k,i,lat)*hox(k,i)/
     |      rmtru_hox(k,i,lat)*oh_h(k,i,lat)*h_hox(k,i,lat)+
     |      gam14(k,i,lat)*o1(k,i)*rmassinv_o1*xnmbar(k,i)/
     |      (.5*(barm(k,i)+barm(k+1,i))))

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

!     call addfld('PDCO2TI' ,' ',' ',pdco2ti,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CO_PROD',' ',' ',co_prod(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CO_LOSS',' ',' ',co_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_co
!-----------------------------------------------------------------------
      subroutine minor_co(tn,o2,ox,n2,he,w,difkk,co,co_nm,hdco,co_out,
     |  conm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use compdat_module,only: wco ! (nlevp1)
      use cons_module,only: rmass_co
      implicit none
!
! Input args:
      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)
     |  ox,       ! oxygen family (mmr)
     |  n2,       ! molecular nitrogen (mmr)
     |  he,      ! helium (mmr)
     |  w,        ! omega (vertical velocity)
     |  difkk,    ! eddy viscosity (from mgw)
     |  co,       ! co (mmr)
     |  co_nm,    ! co at time n-1
     |  hdco      ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  co_out,   ! co output
     |  conm_out  ! co output at time n-1
!
! Local:
      integer :: i,k,ibndbot,ibndtop,nk,nkm1,i0,i1,lat
      real,parameter :: phi_co(3) = (/0.833,1.427,0.852/)
      real,parameter :: alfa_co = 0. ! thermal diffusion coefficient
      real,parameter :: xyco = 1.e-3

      ibndbot = 0 ! ilbc in minor
      ibndtop = 0 ! iubc in minor

      call minor(tn,o2,ox,n2,he,w,difkk,co,co_nm,hdco,co_lbc,
     |  co_ubc,co_loss,co_prod,wco,rmass_co,phi_co,alfa_co,
     |  xyco,ibndbot,ibndtop,co_out,conm_out,
     |  lev0,lev1,lon0,lon1,lat0,lat1,0)

      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!     do lat=lat0,lat1
!       call addfld('CO_OUT' ,' ',' ',co_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('COM_OUT',' ',' ',conm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo

      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 subdomains to boundary conditions:
      allocate(co_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co: error allocating',
     |  ' co_ubc: stat=',i3)") istat
      allocate(co_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co: error allocating',
     |  ' co_lbc: stat=',i3)") istat
      co_ubc = 0. ; co_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(co_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co: error allocating',
     |  ' co_prod: stat=',i3)") istat
      allocate(co_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co: error allocating',
     |  ' co_loss: stat=',i3)") istat
      co_prod = 0. ; co_loss = 0.
!
!     write(6,"('alloc_co: allocated module data')")
!
      end subroutine alloc_co
!-----------------------------------------------------------------------
      end module co_module
