!
      module h2o_module
!
! Advance h2o by one time step.
!
! Boundary conditions, production and loss for H2O are defined
! by comp_h2o, and referenced by minor_h2o. Comp_h2o is called
! from a latitude loop in dynamics. After comp_h2o, dynamics calls
! minor_h2o, 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_h2o (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: h2o_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: h2o_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  h2o_prod,     ! production of h2o (k,i,j)
     |  h2o_loss      ! loss of h2o       (k,i,j)
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_h2o(o1,h2,h2o,h2o2,hox,o1d,op,ch4,cl,phoxic,barm,
     |  xnmbar,lev0,lev1,lon0,lon1,lat)
      use solgar_module,only: xh2olb ! (nlat)
      use cons_module,only: rmass_h2o,rmassinv_h2o2,rmassinv_h2,
     |  rmassinv_ch4,rmassinv_o1
      use qrj_module,only:
     |  pdh2ot    ! total photodissociation of h2o  (F(NJH2OT))
      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: rkm3,rkm30,rkm31,rkm32,rkm33,rkm40,
     |  rkm43,rkm44,rkm45,rk15,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
     |  h2,          ! molecular hydrogen
     |  h2o,         ! water vapor
     |  h2o2,        ! h2o2
     |  hox,         ! hox
     |  o1d,         ! o1(delta)
     |  op,          ! o+
     |  ch4,         ! methane 
     |  cl,          ! cl from comp_solgar
     |  phoxic,      ! from hoxpion
     |  barm,        ! mbar
     |  xnmbar       ! p0e(-z), etc.
!
! Local:
      integer :: i,k,i0,i1,nk,nkm1,ibndbot,ibndtop
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  pdh2oti      ! photodissociation 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 xh2olb.
!
      do i=lon0,lon1 
        h2o_lbc(i,1,lat) = 0.    ! t1
        h2o_lbc(i,2,lat) = 1.    ! t2
        h2o_lbc(i,3,lat) = -xh2olb(lat)*rmass_h2o/barm(lev0,i) ! t3
        h2o_ubc(i,lat) = 0.      ! Zero diffusive flux at top  ! t4
      enddo
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          pdh2oti(k,i) = 0.5*(pdh2ot(k,i,lat)+pdh2ot(k+1,i,lat)) ! s4
!
! Number density production of h2o:
!
          h2o_prod(k,i,lat) = (rkm30(k,i,lat)*oh_h(k,i,lat)**2+
     |      rkm31(k,i,lat)*oh_h(k,i,lat)*ho2_h(k,i,lat)+
     |      rkm40(k,i,lat)*ho2_h(k,i,lat))*(h_hox(k,i,lat)*hox(k,i)/
     |        rmtru_hox(k,i,lat)*xnmbar(k,i))**2+
     |      (rkm33(k,i,lat)*h2(k,i)*rmassinv_h2+
     |        rkm32(k,i,lat)*h2o2(k,i)*rmassinv_h2o2)*xnmbar(k,i)*
     |        oh_h(k,i,lat)*h_hox(k,i,lat)*hox(k,i)/rmtru_hox(k,i,lat)*
     |        xnmbar(k,i)+
     |      (2.*rkm43(k,i,lat)*hox(k,i)/rmtru_hox(k,i,lat)*
     |        oh_h(k,i,lat)*h_hox(k,i,lat)+
     |        rkm44(k,i,lat)*o1(k,i)*rmassinv_o1+
     |        rkm45*o1d(k,i)*rmassinv_o1)*xnmbar(k,i)*ch4(k,i)*
     |        rmassinv_ch4*xnmbar(k,i)+
     |      del3(k,i,lat)*cl(k,i)*xnmbar(k,i)/(.5*(barm(k,i)+
     |        barm(k+1,i)))*ch4(k,i)*rmassinv_ch4*xnmbar(k,i)
!
! Number density loss of h2o:
!
          h2o_loss(k,i,lat) = -(pdh2oti(k,i)+rkm3*o1d(k,i)*rmassinv_o1*
     |      xnmbar(k,i)+rk15*op(k,i)+0.5*phoxic(k,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('PDH2OTI' ,' ',' ',pdh2oti,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2O_PROD',' ',' ',h2o_prod(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('H2O_LOSS',' ',' ',h2o_loss(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_h2o
!-----------------------------------------------------------------------
      subroutine minor_h2o(tn,o2,ox,w,difkk,h2o,h2o_nm,hdh2o,h2o_out,
     |  h2onm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use compdat_module,only: wh2o ! (nlevp1)
      use cons_module,only: rmass_h2o
      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)
     |  h2o,       ! h2o (mmr)
     |  h2o_nm,    ! h2o at time n-1
     |  hdh2o      ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  h2o_out,   ! h2o output
     |  h2onm_out  ! h2o output at time n-1
!
! Local:
      integer :: i,k,ibndbot,ibndtop,nk,nkm1,i0,i1,lat
      real,parameter :: phi_h2o(3) = (/0.817,0.922,0.920/)
      real,parameter :: alfa_h2o = 0. ! thermal diffusion coefficient
      real,parameter :: xyh2o = 1.e-4

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

      call minor(tn,o2,ox,w,difkk,h2o,h2o_nm,hdh2o,h2o_lbc,
     |  h2o_ubc,h2o_loss,h2o_prod,wh2o,rmass_h2o,phi_h2o,alfa_h2o,
     |  xyh2o,ibndbot,ibndtop,h2o_out,h2onm_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('H2O_OUT' ,' ',' ',h2o_out  (:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('H2OM_OUT',' ',' ',h2onm_out(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!     enddo

      end subroutine minor_h2o
!-----------------------------------------------------------------------
      subroutine alloc_h2o(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(h2o_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2o: error allocating',
     |  ' h2o_ubc: stat=',i3)") istat
      allocate(h2o_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2o: error allocating',
     |  ' h2o_lbc: stat=',i3)") istat
!
! Allocate subdomains to production and loss:
      allocate(h2o_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2o: error allocating',
     |  ' h2o_prod: stat=',i3)") istat
      allocate(h2o_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_h2o: error allocating',
     |  ' h2o_loss: stat=',i3)") istat
!
!     write(6,"('alloc_h2o: allocated module data')")
!
      end subroutine alloc_h2o
!-----------------------------------------------------------------------
      end module h2o_module
