!
      module h2_module
!
! Advance h2 by one time step.
!
! Boundary conditions, production and loss for H2O are defined
! by comp_h2, and referenced by minor_h2. Comp_h2 is called
! from a latitude loop in dynamics. After comp_h2, dynamics calls
! minor_h2, 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_h2 (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: h2_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: h2_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  h2_prod,     ! production of h2 (k,i,j)
     |  h2_loss      ! loss of h2       (k,i,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_h2(o1,h2,h2o,hox,o1d,op,ch4,barm,
     |  xnmbar,lev0,lev1,lon0,lon1,lat)
      use solgar_module,only: xh2lb ! (nlat)
      use cons_module,only: rmass_h2,rmassinv_h2o,
     |  rmassinv_ch4,rmassinv_o1
      use qrj_module,only:
     |  pdh2ol    ! lyman-alpha photodissociation of h2o  (F(NJH2OL))
      use comp_meta_module,only: ! (nlevp1,lon0:lon1,lat0:lat1)
     |  oh_h,     ! oh/h  (ratio1)
     |  ho2_h,    ! ho2/h (ratio2)
     |  h_hox,    ! h/hox (ratio3)
     |  rmtru_hox ! "true" mass of hox
      use chemrates_module,only: rkm4,rkm28,rkm33,rkm38,rkm41,rkm5b,rk14
      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
     |  h2,          ! molecular hydrogen
     |  h2o,         ! water vapor
     |  hox,         ! hox
     |  o1d,         ! o1(delta)
     |  op,          ! o+
     |  ch4,         ! methane 
     |  barm,        ! mbar
     |  xnmbar       ! p0e(-z), etc.
!
! Local:
      integer :: i,k,i0,i1,nk,nkm1,ibndbot,ibndtop
! 5/25/06 btf: changed pdh2oli to pdh2olm to indicate midpoints
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pdh2olm      ! lyman-alpha of h2o at midpoints
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! Value at bottom given by specified number density mixing
! ratio xh2lb.
!
      do i=lon0,lon1 
        h2_lbc(i,1,lat) = 0.    ! t1
        h2_lbc(i,2,lat) = 1.    ! t2
        h2_lbc(i,3,lat) = -xh2lb(lat)*rmass_h2/barm(lev0,i) ! t3
        h2_ubc(i,lat) = 0.      ! Zero diffusive flux at top  ! t4
      enddo
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdh2olm(k,i) = 0.5*(pdh2ol(k,i,lat)+pdh2ol(k+1,i,lat)) ! s4
!
! Number density production of h2:
!
          h2_prod(k,i,lat) = xnmbar(k,i)*0.12*pdh2olm(k,i)*
     |      h2o(k,i)*rmassinv_h2o+
     |      (rkm38*ho2_h(k,i,lat)+
     |       rkm41(k,i,lat)*xnmbar(k,i)/(0.5*(barm(k,i)+barm(k+1,i))))*
     |      (h_hox(k,i,lat)*hox(k,i)/rmtru_hox(k,i,lat)*xnmbar(k,i))**2+
     |      rkm5b*ch4(k,i)*rmassinv_ch4*o1d(k,i)*rmassinv_o1*
     |      xnmbar(k,i)*xnmbar(k,i)

!
! Number density loss of h2:
!
          h2_loss(k,i,lat) = -xnmbar(k,i)*(rkm4*o1d(k,i)*rmassinv_o1+
     |      rkm28(k,i,lat)*o1(k,i)*rmassinv_o1+
     |      rkm33(k,i,lat)*hox(k,i)/rmtru_hox(k,i,lat)*oh_h(k,i,lat)*
     |      h_hox(k,i,lat))-rk14*op(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('PDH2OLM' ,' ',' ',pdh2olm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2_PROD',' ',' ',h2_prod(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2_LOSS',' ',' ',h2_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_h2
!-----------------------------------------------------------------------
      subroutine minor_h2(tn,o2,ox,n2,he,w,difkk,h2,h2_nm,hdh2,h2_out,
     |  h2nm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use compdat_module,only: wh2 ! (nlevp1)
      use cons_module,only: rmass_h2
      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)
     |  h2,        ! h2 (mmr)
     |  h2_nm,     ! h2 at time n-1
     |  hdh2       ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  h2_out,   ! h2 output
     |  h2nm_out  ! h2 output at time n-1
!
! Local:
      integer :: i,k,ibndbot,ibndtop,nk,nkm1,i0,i1,lat
      real,parameter :: phi_h2(3) = (/0.226,0.321,0.282/)
      real,parameter :: alfa_h2 = 0. ! thermal diffusion coefficient
      real,parameter :: xyh2 = 1.e-4

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

      call minor(tn,o2,ox,n2,he,w,difkk,h2,h2_nm,hdh2,h2_lbc,
     |  h2_ubc,h2_loss,h2_prod,wh2,rmass_h2,phi_h2,alfa_h2,
     |  xyh2,ibndbot,ibndtop,h2_out,h2nm_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('H2_OUT' ,' ',' ',h2_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('H2M_OUT',' ',' ',h2nm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo

      end subroutine minor_h2
!-----------------------------------------------------------------------
      subroutine alloc_h2(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(h2_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2: error allocating',
     |  ' h2_ubc: stat=',i3)") istat
      allocate(h2_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2: error allocating',
     |  ' h2_lbc: stat=',i3)") istat
      h2_ubc = 0. ; h2_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(h2_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2: error allocating',
     |  ' h2_prod: stat=',i3)") istat
      allocate(h2_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2: error allocating',
     |  ' h2_loss: stat=',i3)") istat
      h2_prod = 0. ; h2_loss = 0.
!
!     write(6,"('alloc_h2: allocated module data')")
!
      end subroutine alloc_h2
!-----------------------------------------------------------------------
      end module h2_module
