!
      subroutine qinite(tn,o2,o1,no,xnmbari,vo2,vo1,vn2,
     |  lev0,lev1,lon0,lon1,lat)
!
! Calculate night-time ionization:
!
      use lbc,only: b,fb ! b(nlonp4,2,2),fb(nlonp4,2)
      use cons_module,only: rmassinv_o2,rmassinv_ox,rmassinv_n2,
     |  rmassinv_no,check_exp
      use qrj_module,only: qo2p,qop,qn2p,qnp, ! (nlevp1,lon0:lon1,lat0:lat1)
     |  qnop,qnoplya,sflux
      use init_module,only: sfeps
      use input_module,only: f107
      use magfield_module,only: rlatm   ! (nlonp4,nlat)
      use chemrates_module,only: disn2p ! (nlevp1,lon0:lon1,lat0:lat1)
      use fields_module,only: tlbc
      use addfld_module,only: addfld
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,o2,o1,no, ! neutral temperature and species (mmr)
     |  xnmbari,     ! p0*e(-z) at interfaces
     |  vo2,vo1,vn2  ! chapman integrals
!
! Local:
      integer :: k,i,n
      real :: sa(3,3),si(3,3),al(3)
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  tni,o2i,o1i,n2i,   ! fields at interfaces
     |  qn_o2,qn_o1,qn_n2, ! nighttime ionization rates
     |  attlya,            ! Attenuation of lyman-alfa at night
     |  tau,
!
! For cosmic ray calculation:
     |  rho,rho1,rho2,rlat1,rho107,f107fac,rlatrho1,rlatrho2,rhofinal
!
! External:
      real,external :: expo
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
#ifdef VT
!     code = 128 ; state = 'qinite' ; activity='ModelCode'
      call vtbegin(128,ier)
#endif
!
! First dimension of 3 is for o2,o,n2, second dim is wavelength:
!               o2       o1        n2
      sa(:,1)=(/1.6E-18, 0.,       0.       /)
      sa(:,2)=(/22.E-18,10.24E-18, 23.11E-18/)
      sa(:,3)=(/16.E-18, 8.40E-18, 11.61E-18/)
      si(:,1)=(/1.0E-18, 0.,       0.       /)
      si(:,2)=(/22.E-18,10.24E-18, 23.11E-18/)
      si(:,3)=(/16.E-18, 8.40E-18, 11.61E-18/)
      al=(/1.5E7,1.5E6,1.5E6/)
!
! o2, o1, n2, tni at interfaces:
      do i=lon0,lon1
        do k=lev0,lev1-1
          o2i(k+1,i) = 0.5*(o2(k,i)+o2(k+1,i))
          o1i(k+1,i) = 0.5*(o1(k,i)+o1(k+1,i))
        enddo
!
! Bottom boundary:
        o2i(1,i) = .5*((b(i,1,1)+1.)*o2(1,i)+ b(i,1,2)*o1(1,i)+
     |    fb(i,1,lat))
        o1i(1,i) = .5* (b(i,2,1)*o2(1,i)+(b(i,2,2)+1.)*o1(1,i)+
     |    fb(i,2,lat))
!
! n2 at interfaces:
        n2i(:,i) = 1.-o2i(:,i)-o1i(:,i)
!
! tn at interfaces:
        do k=lev0+1,lev1
          tni(k,i) = .5*(tn(k,i)+tn(k-1,i))
        enddo ! k=lev0+1,lev1
!       tni(1,i) = tn(lev1,i) ! bottom boundary stored in top slot
        tni(lev0,i) = tlbc(i,lat) ! 5/4/06 btf: lower boundary is in tlbc
        tni(lev1,i) = tn(lev1-1,i)
      enddo ! i=lon0,lon1

!     call addfld('o2i',' ',' ',o2i,'ilev',lev0,lev1,
!    |  'lon',lon0,lon1,lat)
!     call addfld('o1i',' ',' ',o1i,'ilev',lev0,lev1,
!    |  'lon',lon0,lon1,lat)
!     call addfld('n2i',' ',' ',n2i,'ilev',lev0,lev1,
!    |  'lon',lon0,lon1,lat)
!     call addfld('tni',' ',' ',tni,'ilev',lev0,lev1,
!    |  'lon',lon0,lon1,lat)
!     call addfld('xnmbari',' ',' ',xnmbari(:,lon0:lon1),
!    |  'ilev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Init: 
      qn_o2 = 0. ! s8
      qn_o1 = 0. ! s7
      qn_n2 = 0. ! s6
!
! Summation over wavelength (1st dim of sa is species, 2nd is wavelength):
      do n=1,3
        tau = 0.
        do i=lon0,lon1 
          do k=lev0,lev1
            tau(k,i) = tau(k,i)+sa(1,n)*vo2(k,i)
            tau(k,i) = tau(k,i)+sa(2,n)*vo1(k,i)
            tau(k,i) = tau(k,i)+sa(3,n)*vn2(k,i)
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1 
        if (.not.check_exp) then
          tau(:,:) = exp(-tau(:,:))
        else
          tau(:,:) = expo(-tau(:,:))
        endif

!       if (n==1)
!    |    call addfld('tau1',' ',' ',tau,
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       if (n==2)
!    |    call addfld('tau2',' ',' ',tau,
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)
!       if (n==3)
!    |    call addfld('tau3',' ',' ',tau,
!    |      'lev',lev0,lev1,'lon',lon0,lon1,lat)

        do i=lon0,lon1 
          do k=lev0,lev1
            qn_o2(k,i) = qn_o2(k,i)+al(n)*si(1,n)*o2i(k,i)*tau(k,i)*
     |        rmassinv_o2
            qn_o1(k,i) = qn_o1(k,i)+al(n)*si(2,n)*o1i(k,i)*tau(k,i)*
     |        rmassinv_ox
            qn_n2(k,i) = qn_n2(k,i)+al(n)*si(3,n)*n2i(k,i)*tau(k,i)*
     |        rmassinv_n2
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1 
      enddo ! n=1,3

!     call addfld('qn_o2',' ',' ',qn_o2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qn_o1',' ',' ',qn_o1,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qn_n2',' ',' ',qn_n2,
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! Add to total ionization rates:
      do i=lon0,lon1 
        do k=lev0,lev1
          qo2p(k,i,lat) = qo2p(k,i,lat)+0.67*qn_o2(k,i)*xnmbari(k,i)
          qop (k,i,lat) = qop (k,i,lat)+(0.33*qn_o2(k,i)+qn_o1(k,i))*
     |      xnmbari(k,i)
          qn2p(k,i,lat) = qn2p(k,i,lat)+0.86*qn_n2(k,i)*xnmbari(k,i)
          qnp (k,i,lat) = qnp (k,i,lat)+0.14*qn_n2(k,i)*xnmbari(k,i)

        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1 
!
! qnoplya: no+ nighttime ionization by lyman-alpha (XJNOPN):
! (this was moved here from qrj)
      if (.not.check_exp) then
        do i=lon0,lon1 
          do k=lev0,lev1
            attlya(k,i) = (0.68431  *exp(-8.22114e-21*vo2(k,i))+
     |                     0.229841 *exp(-1.77556e-20*vo2(k,i))+
     |                     0.0865412*exp(-8.22112e-21*vo2(k,i)))*sfeps
            qnoplya(k,i,lat) = 2.02e-18*sflux(12)/100.*attlya(k,i)
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1 
      else
        do i=lon0,lon1 
          do k=lev0,lev1
            attlya(k,i) = (0.68431  *expo(-8.22114e-21*vo2(k,i))+
     |                     0.229841 *expo(-1.77556e-20*vo2(k,i))+
     |                     0.0865412*expo(-8.22112e-21*vo2(k,i)))*sfeps
            qnoplya(k,i,lat) = 2.02e-18*sflux(12)/100.*attlya(k,i)
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1 
      endif
 
!     call addfld('qo2p',' ',' ',qo2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qop' ,' ',' ',qop (:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qn2p',' ',' ',qn2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qnp' ,' ',' ',qnp (:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qnoplya',' ',' ',qnoplya(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

!
! Calculate NO ionization and add to NO+
      do i=lon0,lon1 
!
! Bottom boundary:
        qnop(1,i,lat) = qnop(1,i,lat)+qnoplya(1,i,lat)*no(1,i)*
     |    xnmbari(1,i)*rmassinv_no
!
! Levels 2 to top:
        do k=lev0+1,lev1
          qnop(k,i,lat) = qnop(k,i,lat)+qnoplya(k,i,lat)*
     |      0.5*(no(k,i)+no(k-1,i))*xnmbari(k,i)*rmassinv_no
        enddo ! k=lev0+1,lev1
      enddo ! i=lon0,lon1 

!     call addfld('qnop',' ',' ',qnop(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
! Calculate ionization rate from galactic cosmic rays and add to
! total ionization rate:
      do i=lon0,lon1 
        do k=lev0,lev1
          rho(k,i) = (o2i(k,i)*rmassinv_o2+o1i(k,i)*rmassinv_ox+ ! S15
     |                n2i(k,i)*rmassinv_n2)*xnmbari(k,i)
          rho1(k,i) = 1.440e-17*rho(k,i)                         ! S14
          rho2(k,i) = 1.932e-17*rho(k,i)                         ! S13
          rho107(k,i) = rho2(k,i)+(rho1(k,i)-rho2(k,i))/135.*(f107-65.) ! S12
          f107fac(k,i) = 2.84e-17-6.7407e-20*(f107-65.)          ! S14
          rlat1(k,i) = 0.6+0.8*abs(cos(rlatm(i,lat)))            ! S13
          rlatrho1(k,i) =                                        ! S11
     |      (1.74e-18+f107fac(k,i)*(abs(sin(rlatm(i,lat))))**4)*
     |      (3.E+17)**(1.-rlat1(k,i))*rho(k,i)**rlat1(k,i)
          rlatrho2(k,i) =                                        ! S10
     |      (1.74e-18+f107fac(k,i)*(abs(sin(rlatm(i,lat))))**4)*rho(k,i)
          rhofinal(k,i) = rlatrho1(k,i)                          ! S9
          if (rho(k,i) < 3.e+17) rhofinal(k,i) = rlatrho2(k,i)
          if (abs(rlatm(i,lat)) >= 0.925) rhofinal(k,i) = rho107(k,i)
          qo2p(k,i,lat) = qo2p(k,i,lat)+0.154*rhofinal(k,i)
          qop (k,i,lat) = qop (k,i,lat)+0.076*rhofinal(k,i)
          qn2p(k,i,lat) = qn2p(k,i,lat)+0.585*rhofinal(k,i)
!
! disn2p(nlevp1,lon0:lon1,lat0:lat1)
          disn2p(k,i,lat) = disn2p(k,i,lat)+0.585*rhofinal(k,i)
          qnp(k,i,lat) = qnp(k,i,lat)+0.185*rhofinal(k,i)
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1 

!     call addfld('qo2p',' ',' ',qo2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qop' ,' ',' ',qop (:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qn2p',' ',' ',qn2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('disn2p',' ',' ',disn2p(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('qnp',' ',' ',qnp(:,:,lat),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

#ifdef VT
!     code = 128 ; state = 'qinite' ; activity='ModelCode'
      call vtend(128,ier)
#endif
      end subroutine qinite