!
      module ch4_module
!
! Advance ch4 by one time step.
!
! Boundary conditions, production and loss for CH4 are defined
! by comp_ch4, and referenced by minor_ch4. Comp_ch4 is called
! from a latitude loop in dynamics. After comp_ch4, dynamics calls
! minor_ch4, 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_ch4 (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: ch4_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: ch4_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  ch4_prod,     ! production of ch4 (k,i,j)
     |  ch4_loss      ! loss of ch4       (k,i,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_ch4(o1,ch4,hox,o1d,cl,barm,xnmbar,
     |  lev0,lev1,lon0,lon1,lat)
      use solgar_module,only: xch4lb ! (nlat)
      use cons_module,only: rmass_ch4,rmassinv_o1
      use qrj_module,only:
     |  pdch4t  ! total photodissociation of ch4  (F(NJCH4T))
      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: rkm43,rkm44,rkm45,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
     |  ch4,         ! methane from previous step
     |  hox,         ! hox
     |  o1d,         ! o1(delta)
     |  cl,          ! cl from comp_solgar
     |  barm,        ! mbar
     |  xnmbar       ! p0e(-z), etc.
!
! Local:
      integer :: i,k,i0,i1,ibndbot,ibndtop
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pdch4ti,  ! photodissociation of ch4 at midpoints
     |  src_oh    ! OH source
!
      i0=lon0 ; i1=lon1
!
      do i=lon0,lon1 
        ch4_lbc(i,1,lat) = 0.    ! t1
        ch4_lbc(i,2,lat) = 1.    ! t2
        ch4_lbc(i,3,lat) = -xch4lb(lat)*rmass_ch4/barm(lev0,i)
        ch4_ubc(i,lat)   = 0.
      enddo
!     write(6,"('comp_ch4: lat=',i3,' ch4_lbc(3)=',/,(6e12.4))") 
!    |  lat,ch4_lbc(:,3,lat)
!
! Photodiss at midpoints:
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdch4ti(k,i) = 0.5*(pdch4t(k,i,lat)+pdch4t(k+1,i,lat))   ! s4
          src_oh(k,i) = hox(k,i)/rmtru_hox(k,i,lat)*oh_h(k,i,lat)* ! s15
     |      h_hox(k,i,lat)*xnmbar(k,i)
!
! Number density production of ch4:
!
          ch4_prod(k,i,lat) = 0.                                 ! s2
!
! Number density loss of ch4:
!
          ch4_loss(k,i,lat) = -(rkm43(k,i,lat)*src_oh(k,i)+      ! s1
     |      (rkm45*o1d(k,i)*rmassinv_o1+rkm44(k,i,lat)*o1(k,i)*
     |      rmassinv_o1)*xnmbar(k,i)+pdch4ti(k,i)+del3(k,i,lat)*
     |      cl(k,i)*xnmbar(k,i)/(0.5*(barm(k,i)+barm(k+1,i))))

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

!     call addfld('PDCH4TI' ,' ',' ',pdch4ti,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('SRC_OH'  ,' ',' ',src_oh ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CH4_LOSS',' ',' ',ch4_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_ch4
!-----------------------------------------------------------------------
      subroutine minor_ch4(tn,o2,ox,w,difkk,ch4,ch4_nm,hdch4,ch4_out,
     |  ch4nm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use compdat_module,only: wch4 ! (nlevp1)
      use cons_module,only: rmass_ch4
      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)
     |  ch4,     ! ch4 (mmr)
     |  ch4_nm,  ! ch4 at time n-1
     |  hdch4    ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  ch4_out,    ! ch4 output
     |  ch4nm_out   ! ch4 output at time n-1
!
! Local:
      integer :: i,k,ibndbot,ibndtop,i0,i1,lat
      real,parameter :: phi_ch4(3) = (/0.921,0.846,1.077/)
      real,parameter :: alfa_ch4 = 0. ! thermal diffusion coefficient
      real,parameter :: xych4 = 1.e-10

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

      call minor(tn,o2,ox,w,difkk,ch4,ch4_nm,hdch4,ch4_lbc,
     |  ch4_ubc,ch4_loss,ch4_prod,wch4,rmass_ch4,phi_ch4,alfa_ch4,
     |  xych4,ibndbot,ibndtop,ch4_out,ch4nm_out,
     |  lev0,lev1,lon0,lon1,lat0,lat1,0)
!
      i0=lon0 ; i1=lon1
!     do lat=lat0,lat1
!       call addfld('CH4_MNR' ,' ',' ',ch4_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('CH4M_MNR',' ',' ',ch4nm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo

      end subroutine minor_ch4
!-----------------------------------------------------------------------
      subroutine alloc_ch4(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(ch4_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_ch4: error allocating',
     |  ' ch4_ubc: stat=',i3)") istat
      allocate(ch4_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_ch4: error allocating',
     |  ' ch4_lbc: stat=',i3)") istat
      ch4_ubc = 0. ; ch4_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(ch4_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_ch4: error allocating',
     |  ' ch4_prod: stat=',i3)") istat
      allocate(ch4_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_ch4: error allocating',
     |  ' ch4_loss: stat=',i3)") istat
      ch4_prod = 0. ; ch4_loss = 0.
!
!     write(6,"('alloc_ch4: allocated module data')")
!
      end subroutine alloc_ch4
!-----------------------------------------------------------------------
      end module ch4_module
