!
      module co2_module
!
! Advance co2 by one time step.
!
! Boundary conditions, production and loss for CO2 are defined
! by comp_co2, and referenced by minor_co2. Comp_co2 is called
! from a latitude loop in dynamics. After comp_co2, dynamics calls
! minor_co2, which passes this module data to sub minor. Sub
! minor contains 3d mpi calls and its own latitude loops. 
!
      use params_module,only: nlevp1
      implicit none
!
! Boundary conditions and production and loss terms are allocated
! subdomains by sub alloc_co2 (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: co2_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: co2_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  co2_prod,     ! production of co2 (k,i,j)
     |  co2_loss      ! loss of co2       (k,i,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_co2(o1,co,hox,o1d,op,barm,xnmbar,
     |  lev0,lev1,lon0,lon1,lat)
      use solgar_module,only: xco2lb ! (nlat)
      use cons_module,only: rmass_co2,rmassinv_o1,rmassinv_co
      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,rkm42,co2mix,gam14,gam15
      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
     |  hox,         ! hox
     |  o1d,         ! o1(delta)
     |  op,          ! o+
     |  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 
        co2_lbc(i,1,lat) = 0.    ! t1
        co2_lbc(i,2,lat) = 1.    ! t2
        co2_lbc(i,3,lat) = -co2mix*rmass_co2/barm(lev0,i)
        co2_ubc(i,lat)   = 0.
      enddo
!     write(6,"('comp_co2: lat=',i3,' co2_lbc(3)=',/,(6e12.4))") 
!    |  lat,co2_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 co2:
!
          co2_prod(k,i,lat) = xnmbar(k,i)**2*(rkm42(k,i,lat)*     ! s2
     |      co(k,i)*rmassinv_co*hox(k,i)/rmtru_hox(k,i,lat)*
     |      oh_h(k,i,lat)*h_hox(k,i,lat)+gam14(k,i,lat)*co(k,i)*
     |      rmassinv_co*o1(k,i)*rmassinv_o1*xnmbar(k,i)/
     |      (0.5*(barm(k,i)+barm(k+1,i))))
!
! Number density loss of co2:
!
          co2_loss(k,i,lat) = -(rk13*op(k,i)+pdco2ti(k,i)+gam15*   ! s1
     |      o1d(k,i)*rmassinv_o1*xnmbar(k,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('PDCO2TI' ,' ',' ',pdco2ti,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CO2_PROD',' ',' ',co2_prod(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CO2_LOSS',' ',' ',co2_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_co2
!-----------------------------------------------------------------------
      subroutine minor_co2(tn,o2,ox,w,difkk,co2,co2_nm,hdco2,co2_out,
     |  co2nm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use compdat_module,only: wco2 ! (nlevp1)
      use cons_module,only: rmass_co2
      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)
     |  w,       ! omega (vertical velocity)
     |  difkk,   ! eddy viscosity (from mgw)
     |  co2,     ! co2 (mmr)
     |  co2_nm,  ! co2 at time n-1
     |  hdco2    ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  co2_out,    ! co2 output
     |  co2nm_out   ! co2 output at time n-1
!
! Local:
      integer :: i,k,ibndbot,ibndtop,nk,nkm1,i0,i1,lat
      real,parameter :: phi_co2(3) = (/1.199,3.91,1.201/)
      real,parameter :: alfa_co2 = 0. ! thermal diffusion coefficient
      real,parameter :: xyco2 = 1.e-4

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

!     subroutine minor(tn,o2,ox,w,difkk,fin,fnm,fhd,flbc,fubc,
!    |  sloss,sprod,small,rmass,phix,alfax,xy,ilbc,iubc,fout,fout_nm,
!    |  lev0,lev1,lon0,lon1,lat0,lat1,idebug)

      call minor(tn,o2,ox,w,difkk,co2,co2_nm,hdco2,co2_lbc,
     |  co2_ubc,co2_loss,co2_prod,wco2,rmass_co2,phi_co2,alfa_co2,
     |  xyco2,ibndbot,ibndtop,co2_out,co2nm_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('CO2_OUT' ,' ',' ',co2_out  (:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CO2M_OUT',' ',' ',co2nm_out(:,lon0:lon1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo

      end subroutine minor_co2
!-----------------------------------------------------------------------
      subroutine alloc_co2(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(co2_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co2: error allocating',
     |  ' co2_ubc: stat=',i3)") istat
      allocate(co2_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co2: error allocating',
     |  ' co2_lbc: stat=',i3)") istat
      co2_ubc = 0. ; co2_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(co2_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co2: error allocating',
     |  ' co2_prod: stat=',i3)") istat
      allocate(co2_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_co2: error allocating',
     |  ' co2_loss: stat=',i3)") istat
      co2_prod = 0. ; co2_loss = 0.
!
!     write(6,"('alloc_co2: allocated module data')")
!
      end subroutine alloc_co2
!-----------------------------------------------------------------------
      end module co2_module
