!
      module n4s_module
!
! Boundary conditions, production and loss for N4S are defined
! by comp_n4s, and referenced by minor_n4s. Comp_n4s is called
! from a latitude loop in dynamics. After comp_n4s, dynamics calls
! minor_n4s, 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 compdat_module,only: wn4s ! (nlevp1)
      use addfld_module,only: addfld
      implicit none
!
! Boundary conditions and production and loss terms are allocated
! subdomains by sub alloc_n4s (called from allocdata.F).
!
      real,allocatable,dimension(:,:)   :: n4s_ubc ! upper boundary (i,j)
      real,allocatable,dimension(:,:,:) :: n4s_lbc ! lower boundary (i,3,j)
      real,allocatable,dimension(:,:,:) ::
     |  n4s_prod,     ! production of n4s (k,i,j)
     |  n4s_loss      ! loss of n4s       (k,i,j)
      real :: phi_n4s(3) = (/0.651, 0.731, 0.741/)
      real,parameter :: alfa_n4s = 0. ! thermal diffusion coefficient
!
      contains
!-----------------------------------------------------------------------
      subroutine comp_n4s(tn,o2,ox,o1,no,no2,oh,hox,ne,n2d,barm,xnmbar,
     |  o2p,op,n2p,nplus,nop,lev0,lev1,lon0,lon1,lat)
!
! Advance n4s.
!
      use qrj_module,only: 
     |  pdn2,   ! f(nqtef)
     |  pdnosrb ! xjno
      use cons_module,only: p0,expz,boltz,expzmid_inv,brn2d,
     |  rmassinv_n2d,rmassinv_o1,rmassinv_no,rmassinv_n2,rmassinv_o2,
     |  rmassinv_no2,rmass_n4s
      use chemrates_module,only: beta1,beta3,beta4,beta5,beta7,beta8,
     |  beta13,rk2,rk4,rk6,rk8,ra1,ra3
      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 fields_module,only: tlbc ! lower boundary interface level of TN
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn       ,   ! tn (deg K) 
     |  o2       ,   ! o2 (mmr)
     |  ox       ,   ! ox (mmr) (not used?)
     |  o1       ,   ! o1 (mmr)
     |  no       ,   ! no (mmr)
     |  no2      ,   ! no2(mmr)
     |  oh       ,   ! oh (mmr)
     |  hox      ,   ! hox(mmr)
     |  ne       ,   ! ne (cm3)
     |  n2d      ,   ! updated n2d from comp_meta
     |  barm     ,   ! mbar
     |  xnmbar   ,   ! n*mbar
     |  o2p      ,   ! o2+
     |  op       ,   ! o+
     |  n2p      ,   ! n2+
     |  nplus    ,   ! n+
     |  nop          ! no+
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1
      real :: xnmbarlb,n2
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! Check inputs:
!     call addfld('TN_N4S'  ,' ',' ',tn  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O2_N4S'  ,' ',' ',o2  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OX_N4S'  ,' ',' ',ox  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O1_N4S'  ,' ',' ',o1  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NO_N4S'  ,' ',' ',no  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NO2_N4S' ,' ',' ',no2 (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OH_N4S'  ,' ',' ',oh  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('HOX_N4S' ,' ',' ',hox (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NE_N4S'  ,' ',' ',ne  (:,i0:i1), 
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)           ! Ne is at interfaces
!     call addfld('N2D_N4S' ,' ',' ',n2d (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('MBAR_N4S',' ',' ',barm(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O2P_N4S' ,' ',' ',o2p (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OP_N4S'  ,' ',' ',op  (lev0:lev1-1,i0:i1),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('N2P_N4S' ,' ',' ',n2p (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NOP_N4S' ,' ',' ',nop (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NP_N4S'  ,' ',' ',nplus(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Lower boundary:
! n4s_lbc(:,1)=A, n4s_lbc(:,2)=B, n4s_lbc(:,3)=C define lower boundary
!   condition where A*DPSX/DZ + B*PSX + C = 0.
! Note tn(lev1,i) assumes tn bottom boundary is stored in top slot.
! 5/4/06 btf: bottom boundary of TN is in tlbc:
!
      do i=lon0,lon1 
!       xnmbarlb = p0*expz(lev0)*expzmid_inv*barm(lev0,i)/  ! t5
!    |    (boltz*tn(lev1,i))
        xnmbarlb = p0*expz(lev0)*expzmid_inv*barm(lev0,i)/
     |    (boltz*tlbc(i,lat))
        n2 = (1.-o2(lev0,i)-o1(lev0,i))
        n4s_lbc(i,1,lat) = 0.
        n4s_lbc(i,2,lat) = 1.
!
        n4s_lbc(i,3,lat) = -(pdn2(lev0,i,lat)*(1.-brn2d)+xnmbarlb*    ! t3
     |    (n2d(lev0,i)*rmassinv_n2d*(xnmbarlb*beta4*o1(lev0,i)*
     |    rmassinv_o1+beta5(lev0,i,lat)*sqrt(ne(lev0,i)+ne(lev0+1,i))+
     |    beta7)+0.5*(pdnosrb(lev0,i,lat)+pdnosrb(lev0+1,i,lat))*
     |    no(lev0,i)*rmassinv_no)+xnmbarlb*(rk2(lev0,i,lat)*op(lev0,i)*
     |    n2*rmassinv_n2+rk6*nplus(lev0,i)*o2(lev0,i)*rmassinv_o2+rk8*
     |    nplus(lev0,i)*o1(lev0,i)*rmassinv_o1)+sqrt(ne(lev0,i)*
     |    ne(lev0+1,i))*(ra1(lev0,i,lat)*nop(lev0,i)*0.15+
     |    ra3(lev0,i,lat)*n2p(lev0,i)*1.1))/(xnmbarlb*(beta1(lev0,i,lat)
     |    *o2(lev0,i)*rmassinv_o2+beta3(lev0,i,lat)*no(lev0,i)*
     |    rmassinv_no+beta8*hox(lev0,i)/rmtru_hox(lev0,i,lat)*
     |    oh_h(lev0,i,lat)*h_hox(lev0,i,lat)+beta13*no2(lev0,i)*
     |    rmassinv_no2)-rk4*no2(lev0,i))*rmass_n4s/xnmbarlb
!
! Zero diffusive flux at top:
        n4s_ubc(i,lat) = 0.
      enddo ! i=lon0,lon1 

!     write(6,"('comp_n4s: lat=',i3,' n4s_lbc(i,3,lat)=',/,(6e12.4))")
!    |  lat,n4s_lbc(:,3,lat)
!
! Sources:
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! Number density production:
! (Small diamond diffs with tgcm24 at lbc)
          n2 = (1.-o2(k,i)-o1(k,i))
          n4s_prod(k,i,lat) = 0.5*(pdn2(k,i,lat)+pdn2(k+1,i,lat))*      ! s2
     |      (1.-brn2d)+xnmbar(k,i)*(n2d(k,i)*rmassinv_n2d*(xnmbar(k,i)*
     |      beta4*o1(k,i)*rmassinv_o1+beta5(k,i,lat)*sqrt(ne(k,i)+
     |      ne(k+1,i))+beta7)+0.5*(pdnosrb(k,i,lat)+pdnosrb(k+1,i,lat))*
     |      no(k,i)*rmassinv_no)+xnmbar(k,i)*(rk2(k,i,lat)*op(k,i)*n2*
     |      rmassinv_n2+rk6*nplus(k,i)*o2(k,i)*rmassinv_o2+rk8*
     |      nplus(k,i)*o1(k,i)*rmassinv_o1)+sqrt(ne(k,i)*ne(k+1,i))*
     |      (ra1(k,i,lat)*nop(k,i)*0.15+ra3(k,i,lat)*n2p(k,i)*1.1)
!
! Number density loss:
          n4s_loss(k,i,lat) = -xnmbar(k,i)*(beta1(k,i,lat)*o2(k,i)*     ! s1
     |      rmassinv_o2+beta3(k,i,lat)*no(k,i)*rmassinv_no+beta8*
     |      hox(k,i)/rmtru_hox(k,i,lat)*oh_h(k,i,lat)*h_hox(k,i,lat)+
     |      beta13*no2(k,i)*rmassinv_no2)-rk4*o2p(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('N4S_PROD'  ,' ',' ',n4s_prod(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('N4S_LOSS'  ,' ',' ',n4s_loss(:,i0:i1,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      end subroutine comp_n4s
!-----------------------------------------------------------------------
      subroutine minor_n4s(tn,o2,ox,w,difkk,n4s,n4s_nm,hdn4s,n4s_out,
     |  n4snm_out,lev0,lev1,lon0,lon1,lat0,lat1)
      use cons_module,only: rmass_n4s
!
! Call minor to complete n4s. This is called from dynamics, once
! per time step after comp_n4s (which is called from a latitude loop).
! Inputs are at 3d subdomains. 
!
      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)
     |  n4s,     ! n4s (mmr)
     |  n4s_nm,  ! n4s at time n-1
     |  hdn4s    ! horizontal diffusion
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  n4s_out,    ! n4s output
     |  n4snm_out   ! new n4s at time n-1
!
! Local:
      integer :: ibndbot,ibndtop,nk,nkm1,i0,i1
      real,parameter :: phin4s(3)=(/0.651,0.731,0.741/)
      real,parameter :: xyn4s = 1.e-6

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

      call minor(tn,o2,ox,w,difkk,n4s,n4s_nm,hdn4s,n4s_lbc,
     |  n4s_ubc,n4s_loss,n4s_prod,wn4s,rmass_n4s,phin4s,alfa_n4s,
     |  xyn4s,ibndbot,ibndtop,n4s_out,n4snm_out,
     |  lev0,lev1,lon0,lon1,lat0,lat1,0)

      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1

      end subroutine minor_n4s
!-----------------------------------------------------------------------
      subroutine alloc_n4s(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(n4s_ubc(lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_n4s: error allocating',
     |  ' n4s_ubc: stat=',i3)") istat
      allocate(n4s_lbc(lon0:lon1,3,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_n4s: error allocating',
     |  ' n4s_lbc: stat=',i3)") istat
      n4s_ubc = 0. ; n4s_lbc = 0.
!
! Allocate subdomains to production and loss:
      allocate(n4s_prod(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_n4s: error allocating',
     |  ' n4s_prod: stat=',i3)") istat
      allocate(n4s_loss(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_n4s: error allocating',
     |  ' n4s_prod: stat=',i3)") istat
      n4s_prod = 0. ; n4s_loss = 0.

!     write(6,"('alloc_n4s: allocated module data')")
!
      end subroutine alloc_n4s
!-----------------------------------------------------------------------
      endmodule n4s_module