! 
      subroutine settei(tn,o2,o1,ne,te,ti,op,o2p,nop,barm,qji_ti,
     |  te_out,ti_out,lev0,lev1,lon0,lon1,lat)
!
! Calculate electron and ion temperatures.
!
      use params_module,only: dz,nlonp4,spval
      use cons_module,only: pi,rtd,evergs,p0,expz,expzmid,expzmid_inv,
     |  boltz,rmassinv_o2,rmassinv_o1,rmassinv_n2,gask,grav,dipmin,avo
      use input_module,only: f107
      use chapman_module,only: chi   ! solar zenith angle (nlonp4,nlat)
      use magfield_module,only: rlatm,dipmag
      use aurora_module,only: qteaur ! (nlonp4,nlat)
      use qrj_module,only: ! Q is modified, all others are input.
     |  qtotal,! total heating 
     |  qop2p, ! o+(2p)
     |  qop2d, ! o+(2d)
     |  qo2p,  ! o2+ ionization
     |  qop,   ! o+  ionization
     |  qn2p,  ! n2+ ionization
     |  qnp,   ! n+  ionization
     |  qnop   ! no+ ionization
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,      ! neutral temperature (deg K)
     |  o2,      ! molecular oxygen (mmr)
     |  o1,      ! atomic oxygen (mmr)
     |  ne,      ! electron density (cm3)
     |  te,      ! electron temperature (from previous time step)
     |  ti,      ! ion temperature (from previous time step)
     |  op,      ! O+
     |  o2p,     ! O2+
     |  nop,     ! NO+
     |  barm,    ! mean molecular weight
     |  qji_ti   ! joule heating from qjoule_ti (used ui,vi)
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  te_out,  ! output electron temperature (deg K) 
     |  ti_out   ! output ion temperature (deg K)
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: k,i,n,ier
      integer :: nk,nkm1
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  te_int,    ! electron temperature (interfaces)
     |  tn_int,    ! neutral temperature (interfaces)
     |  o2n,       ! O2 number density (midpoints or interfaces)
     |  o1n,       ! O1 number density (midpoints or interfaces)
     |  n2n,       ! N2 number density (midpoints or interfaces)
     |  xnmbar,    ! p0*e(-z)*barm/kT  (minpoints or interfaces)
     |  root_te,   ! sqrt(te)
     |  root_tn,   ! sqrt(tn)
     |  root_ne,   ! sqrt(ne)
     |  tek0,      ! ke/te**2.5 (s15)
     |  h_mid,h_int,
     |  p_coef,    ! coefficient for trisolv     (s1)
     |  q_coef,    ! coefficient for trisolv     (s2)
     |  r_coef,    ! coefficient for trisolv     (s3)
     |  rhs,       ! right-hand-side for trisolv (s4)
     |  qtot,      ! total ionization rate       (s11)
     |  fki,       ! work array
     |  qe,        ! source term                 (s10)
     |  q_eni,     ! heating from electron/neutral and electron/ion collisions
     |  coll_en2v, ! electron/N2vib collision    (s9)
! Lei
     |  qout,      ! electron heating
     |  lossout,   ! electron heating
!
! Cooling rates (heat loss):
     |  loss_en2v, ! electron/N2vib loss term    (s10)
     |  loss_en2,  ! electron/N2 loss
     |  loss_eo2,  ! electron/O2 loss
     |  loss_eo1d, ! electron/O(1d) loss
     |  loss_eo1,  ! electron/O loss
     |  loss_xen,  ! L0*(E,N) (s8)
     |  loss_en,   ! electrons/neutrals loss     (s11)
     |  loss_ei,   ! electron/ion loss           (s10)
     |  loss_in    ! ion/neutral loss            (s9)
      real,parameter ::
     |  fpolar = -3.0e+9,  ! polar te flux
     |  del    = 1.e-6  ,
!
! Correction factors for neutral heating due to L(E,O1D)
     |  alam = 0.0069   ,
     |  ad   = 0.0091   ,
     |  sd   = 2.3e-11
      real :: 
     |  root2,  ! sqrt(2)
     |  f107te  ! solar flux
!
! a,fed,fen,fe,sindipmag have a z dimension only for diagnostic plotting:
      real,dimension(lon0:lon1) ::
     |  a,fed,fen, ! day/night
     |  fe,        ! heat flux at upper boundary
     |  sindipmag  ! sin(dipmag)
!
! For diagnostic plotting (addfsech):
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  a_ki,        ! for diagnostic plotting of a
     |  fed_ki,      ! for diagnostic plotting of fed
     |  fen_ki,      ! for diagnostic plotting of fen
     |  fe_ki,       ! for diagnostic plotting of fe
     |  dipmag_ki,   ! for diagnostic plotting of dipmag
     |  chi_ki,      ! for diagnostic plotting of chi
     |  qteaur_ki,   ! for diagnostic plotting of qteaur
     |  sindipmag_ki ! for diagnostic plotting of sindipmag
!
      real,dimension(lev0:lev1) :: expzp ! e(-z) (midpoints or interfaces)
      real :: bb
!
#ifdef VT
!     code = 126 ; state = 'settei' ; activity='ModelCode'
      call vtbegin(126,ier)
#endif
!
      root2 = sqrt(2.)
      f107te = f107
      if (f107te > 235.) f107te = 235.
      nk = lev1-lev0+1
      nkm1 = nk-1

!     call addfsech('QJI_TI',' ',' ',qji_ti,lon0,lon1,nk,nkm1,lat)
!
      do i=lon0,lon1
        if (abs(rlatm(i,lat))-pi/4.5 >= 0.) then
          a(i) = 1.
        else
          a(i) = .5*(1.+sin(pi*(abs(rlatm(i,lat))-pi/9.)/(pi/4.5)))
        endif 
!
! 8/28/06 btf: mods by Jiuhou Lei:
       if (abs(rlatm(i,lat))-pi/4.5 >= 0.) then
          bb = 5.0
        else
          bb=abs(rlatm(i,lat))*180/pi
          bb = 5.0 -4.3e-3*bb*bb + 6.85
        endif 
!
! Increased heat flux for TE fom protonosphere.
! 8/28/06 btf: mods by Jiuhou:
        fed(i) = -4.0e+7*f107te*a(i)-2.0e+7*f107te   
        fen(i) = fed(i)/bb
!
        fed(i) = fed(i)+qteaur(i,lat)     ! t4
        fen(i) = fen(i)+qteaur(i,lat)     ! t5
        if (chi(i,lat)-.5*pi >= 0.) then  ! chi==t2
          fe(i) = fen(i)                ! t1
        else
          fe(i) = fed(i)
        endif
        if ((chi(i,lat)*rtd-80.)*(chi(i,lat)*rtd-100.)>=0.) then
          fe(i) = fe(i)*evergs
        else
          fe(i) = (.5*(fed(i)+fen(i))+.5*(fed(i)-fen(i))*
     |      cos(pi*(chi(i,lat)*rtd-80.)/20.))*evergs
        endif
!
! Add fpolar if magnetic latitude >= 60 degrees:
        if (abs(rlatm(i,lat))-pi/3.>=0.) fe(i) = fe(i)+fpolar*evergs
!
! For plotting:
        a_ki  (:,i) = a(i)
        chi_ki(:,i) = chi(i,lat)
        qteaur_ki(:,i) = qteaur(i,lat)
        fed_ki(:,i) = fed(i)
        fen_ki(:,i) = fen(i)
        fe_ki (:,i) = fe(i)
      enddo ! i=lon0,lon1

!     call addfsech('MAGLAT',' ',' ',a_ki     ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('CHI'   ,' ',' ',chi_ki   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('QTEAUR',' ',' ',qteaur_ki,lon0,lon1,nk,nkm1,lat)
!     call addfsech('FED'   ,' ',' ',fed_ki   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('FEN'   ,' ',' ',fen_ki   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('FE'    ,' ',' ',fe_ki    ,lon0,lon1,nk,nkm1,lat)

!     write(6,"('settei: lat=',i2,' fed=',/,(6e12.4))") lat,fed
!     write(6,"('settei: lat=',i2,' fen=',/,(6e12.4))") lat,fen
!     write(6,"('settei: lat=',i2,' fe =',/,(6e12.4))") lat,fe

!
! te,o2,o,tn at interfaces: 
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          te_int(k,i) = .5*(te(k,i)+te(k-1,i)) ! s14
          o2n(k,i)    = .5*(o2(k,i)+o2(k-1,i)) ! s13
          o1n(k,i)    = .5*(o1(k,i)+o1(k-1,i)) ! s12
          tn_int(k,i) = .5*(tn(k,i)+tn(k-1,i)) ! s10
        enddo ! k=lev0+1,lev1-2
!
! Bottom:
        te_int(lev0,i) = 1.5*te(lev0,i)-.5*te(lev0+1,i)
        o2n(lev0,i)    = 1.5*o2(lev0,i)-.5*o2(lev0+1,i)
        o1n(lev0,i)    = 1.5*o1(lev0,i)-.5*o1(lev0+1,i)
        tn_int(lev0,i) = 1.5*tn(lev0,i)-.5*tn(lev0+1,i)
!
! Top:
        te_int(lev1,i) = 1.5*te(lev1-1,i)-.5*te(lev1-2,i)
        o2n(lev1,i)    = 1.5*o2(lev1-1,i)-.5*o2(lev1-2,i)
        o1n(lev1,i)    = 1.5*o1(lev1-1,i)-.5*o1(lev1-2,i)
        tn_int(lev1,i) = 1.5*tn(lev1-1,i)-.5*tn(lev1-2,i)
!
! N2:
        do k=lev0,lev1
          n2n(k,i) = 1.-o2n(k,i)-o1n(k,i)
          if (n2n(k,i) < 0.) n2n(k,i) = 0.
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfsech('TE_INT'  ,' ',' ',te_int,lon0,lon1,nk,nkm1,lat)
!     call addfsech('O2_INT'  ,' ',' ',o2n   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('O1_INT'  ,' ',' ',o1n   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('TN_INT'  ,' ',' ',tn_int,lon0,lon1,nk,nkm1,lat)
!     call addfsech('N2_INT'  ,' ',' ',n2n   ,lon0,lon1,nk,nkm1,lat)
!
! Convert o2,o,n2 to number density (interfaces):
      do k=lev0,lev1-1
        expzp(k) = expzmid_inv*expz(k)
      enddo ! k=lev0,lev1
      expzp(lev1) = expzmid*expz(lev1-1)
      do i=lon0,lon1
        do k=lev0,lev1
! xnmbar at interfaces:
          xnmbar(k,i) = p0*expzp(k)*barm(k,i)/(boltz*tn_int(k,i)) ! s8
          o2n(k,i) = xnmbar(k,i)*o2n(k,i)*rmassinv_o2       ! s13
          o1n(k,i) = xnmbar(k,i)*o1n(k,i)*rmassinv_o1       ! s12
          n2n(k,i) = xnmbar(k,i)*n2n(k,i)*rmassinv_n2       ! s11
          root_te(k,i) = sqrt(te_int(k,i))
!
          tek0(k,i) = 7.5e5/
     |      (1.+3.22e4*te_int(k,i)**2/
     |      ne(k,i)*(root_te(k,i)*
     |      (2.82e-17 - 3.41e-21   * te_int (k,i))*n2n(k,i)+
     |      (2.20e-16 + 7.92e-18   * root_te(k,i))*o2n(k,i)+
     |       1.10e-16 * (1.+5.7e-4 * te_int (k,i))*o1n(k,i)))*evergs
!
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfsech('XNMBARI',' ',' ',xnmbar,lon0,lon1,nk,nk,lat) ! s8
!     call addfsech('O2_INT' ,' ',' ',o2n   ,lon0,lon1,nk,nk,lat) ! s13
!     call addfsech('O1_INT' ,' ',' ',o1n   ,lon0,lon1,nk,nk,lat) ! s12
!     call addfsech('N2_INT' ,' ',' ',n2n   ,lon0,lon1,nk,nk,lat) ! s11
!     call addfsech('TEK0'   ,' ',' ',tek0  ,lon0,lon1,nk,nk,lat) ! s15

      do i=lon0,lon1
        do k=lev0,lev1-1
          h_mid(k,i) = gask*tn(k,i)/(.5*(barm(k,i)+barm(k+1,i))*grav) ! s7
        enddo ! k=lev0,lev1-1
        do k=lev0,lev1
          h_int(k,i) = gask*tn_int(k,i)/(barm(k,i)*grav)              ! s6
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfsech('H_MID'   ,' ',' ',h_mid,lon0,lon1,nk,nkm1,lat) ! s7
!     call addfsech('H_INT'   ,' ',' ',h_int,lon0,lon1,nk,nk  ,lat) ! s6


      do i=lon0,lon1
        if (abs(dipmag(i,lat)) >= dipmin) then
          sindipmag(i) = (sin(dipmag(i,lat)))**2 ! t2,s2
        else
          sindipmag(i) = (sin(dipmin))**2
        endif
        if (sindipmag(i) < .10) sindipmag(i) = .10
!
! Start coefficients and rhs for trsolv:
        do k=lev0,lev1-1
          p_coef(k,i) = 2./7.*sindipmag(i)/(h_mid(k,i)*dz**2) ! s1
          r_coef(k,i) = p_coef(k,i)*tek0(k+1,i)/h_int(k+1,i)  ! s3
          p_coef(k,i) = p_coef(k,i)*tek0(k  ,i)/h_int(k  ,i)  ! s1
          q_coef(k,i) = -(p_coef(k,i)+r_coef(k,i))            ! s2
          rhs(k,i) = 0.                                       ! s4
        enddo ! k=lev0,lev1-1
!
! Bottom boundary:
        q_coef(lev0,i) = q_coef(lev0,i)-p_coef(lev0,i)
        rhs(lev0,i) = rhs(lev0,i)-2.*p_coef(lev0,i)*tn_int(lev0,i)**3.5
        p_coef(lev0,i) = 0.
!
! Upper boundary:
        q_coef(lev1-1,i) = q_coef(lev1-1,i)+r_coef(lev1-1,i)
        rhs(lev1-1,i) = rhs(lev1-1,i)+r_coef(lev1-1,i)*dz*3.5*
     |    h_int(lev1,i)*fe(i)/tek0(lev1,i)
        r_coef(lev1-1,i) = 0.
      enddo ! i=lon0,lon1

      do i=lon0,lon1
        dipmag_ki(:,i) = dipmag(i,lat)
        sindipmag_ki(:,i) = sindipmag(i)
      enddo

!     call addfsech('DIPMAG' ,' ',' ',dipmag_ki,lon0,lon1,nk,nkm1,lat)
!     call addfsech('SINDIPM',' ',' ',sindipmag_ki,lon0,lon1,nk,nkm1,
!    |  lat)

!     call addfsech('P_COEF' ,' ',' ',p_coef,lon0,lon1,nk,nkm1,lat) ! s1
!     call addfsech('Q_COEF' ,' ',' ',q_coef,lon0,lon1,nk,nkm1,lat) ! s2
!     call addfsech('R_COEF' ,' ',' ',r_coef,lon0,lon1,nk,nkm1,lat) ! s3
!     call addfsech('RHS0'   ,' ',' ',rhs   ,lon0,lon1,nk,nkm1,lat) ! s4

!
! qtot = total ionization rate = sum(Qxx) = 
! (QO2+) + (QO+) + (QN2+) + (QNO+) + (QN+) + (QO+(2D)) + (QO+(2P))
!

!     call addfsech('QO2P' ,' ',' ',qo2p (:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QOP'  ,' ',' ',qop  (:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QN2P' ,' ',' ',qn2p (:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QNOP' ,' ',' ',qnop (:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QNP'  ,' ',' ',qnp  (:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QOP2D',' ',' ',qop2d(:,:,lat),lon0,lon1,nk,nk,lat)
!     call addfsech('QOP2P',' ',' ',qop2p(:,:,lat),lon0,lon1,nk,nk,lat)

!
! Note re comparison with tgcm15: Because the check for qtot < 1.e-20
! is inside the reduction operation loop, tiegcm1 must sum separately
! for each ion species, as in tgcm15. If tiegcm1 sums all species at
! each grid point (qtot(k,i)=qo2p(k,i)+qop(k,i)+...), there are diffs
! at the bottom boundary (which is where qtot < 1.e-20). By summing
! each species separately, as below, there are no diffs. Ions are also
! summed in qjion.F, but there is no check for < 1.e-20 there, so summing
! all species at each grid point works fine (see qtot in qjion.F).
!
      qtot = 0. ! whole array init
      do n=1,7
        select case (n)
          case (1) 
            fki(:,:) = qo2p(:,:,lat)
          case (2) 
            fki(:,:) = qop(:,:,lat)
          case (3) 
            fki(:,:) = qn2p(:,:,lat)
          case (4) 
            fki(:,:) = qnop(:,:,lat)
          case (5) 
            fki(:,:) = qnp(:,:,lat)
          case (6) 
            fki(:,:) = qop2d(:,:,lat)
          case (7) 
            fki(:,:) = qop2p(:,:,lat)
        end select
        do i=lon0,lon1
          do k=lev0,lev1
            qtot(k,i) = qtot(k,i)+fki(k,i)
            if (qtot(k,i) < 1.e-20) qtot(k,i) = 1.e-20
          enddo
        enddo
      enddo ! n=1,7

!     call addfsech('QTOT_SUM',' ',' ',qtot,lon0,lon1,nk,nk,lat) ! s11

      do i=lon0,lon1
        do k=lev0,lev1-1
          qtot(k,i) = sqrt(qtot(k,i)*qtot(k+1,i))
        enddo
        qtot(lev1,i) = 0.
      enddo ! i=lon0,lon1

!     call addfsech('QTOT',' ',' ',qtot   ,lon0,lon1,nk,nk,lat) ! s11

!
! Set Ne (midpoints "(K+1/2)"):
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          root_ne(k,i) = ne(k,i)*ne(k+1,i)
          if (root_ne(k,i) < 1.e4) root_ne(k,i) = 1.e4
          root_ne(k,i) = sqrt(root_ne(k,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfsech('ROOT_NE',' ',' ',root_ne,lon0,lon1,nk,nkm1,lat) ! s15
!
! Set up o2,o,n2 number densities at midpoints:
!
      do i=lon0,lon1 ! (DO 20)
        do k=lev0,lev1-1
          xnmbar(k,i) = p0*expz(k)*.5*(barm(k,i)+barm(k+1,i))/
     |      (boltz*tn(k,i))
          o2n(k,i) = xnmbar(k,i)*o2(k,i)*rmassinv_o2  ! s14
          o1n(k,i) = xnmbar(k,i)*o1(k,i)*rmassinv_o1  ! s13
          n2n(k,i) = 1.-o2(k,i)-o1(k,i)
          if (n2n(k,i) < 0.) n2n(k,i) = 0.
          n2n(k,i) = xnmbar(k,i)*n2n(k,i)*rmassinv_n2 ! s12
!
! Calculate source term qe (s10)
! Comment from earlier version (maybe the *1.0 below was once *2.0):
!   "Correction facor of 2 increase in TE heating rate"
!
          qe(k,i) = alog(root_ne(k,i)/(o2n(k,i)+n2n(k,i)+0.1*o1n(k,i)))
          qe(k,i) = exp(-((((0.001996*qe(k,i)+0.08034)*qe(k,i)+1.166)*
     |              qe(k,i)+6.941)*qe(k,i)+12.75))*1.0
!
! Subtract qe from right-hand-side:
          rhs(k,i) = rhs(k,i)-qe(k,i)*qtot(k,i)*evergs
          root_te(k,i) = sqrt(te(k,i))
! by Lei
          qout(k,i) = qe(k,i)*qtot(k,i)*evergs

!
! Electron/N2 collision A(E,N2,VIB) (s9):
!
        if (te(k,i) >= 1000.) then
          coll_en2v(k,i) = 2.e-7*exp(-4605.2/te(k,i))
        else
          coll_en2v(k,i) = 5.71e-8*exp(-3352.6/te(k,i))
        endif
        if (te(k,i) > 2000.)
     |    coll_en2v(k,i) = 2.53e-6*root_te(k,i)*exp(-17620./te(k,i))
!
! Loss due to electron/n2 collision L0(E,N2,VIB)/(NE*N(N2)) (s10)
!
          loss_en2v(k,i) = 3200.*(1./te(k,i)-1./tn(k,i))
          loss_en2v(k,i) = sign(abs(loss_en2v(k,i))+del,loss_en2v(k,i))
          loss_en2v(k,i) = -3200./(te(k,i)*tn(k,i))*
     |      (1.-exp(loss_en2v(k,i)))/loss_en2v(k,i)
          loss_en2v(k,i) = 1.3e-4*loss_en2v(k,i)*coll_en2v(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1 (end DO 20)

!     call addfsech('XNMBARM',' ',' ',xnmbar,lon0,lon1,nk,nkm1,lat) ! s6
!     call addfsech('O2N'    ,' ',' ',o2n   ,lon0,lon1,nk,nkm1,lat) ! s14
!     call addfsech('O1N'    ,' ',' ',o1n   ,lon0,lon1,nk,nkm1,lat) ! s13
!     call addfsech('N2N'    ,' ',' ',n2n   ,lon0,lon1,nk,nkm1,lat) ! s12

!     call addfsech('C_EN2V'  ,' ',' ',coll_en2v,lon0,lon1,nk,nkm1,lat) ! s9
!     call addfsech('L_EN2V'  ,' ',' ',loss_en2v,lon0,lon1,nk,nkm1,lat) ! s10

!
! Calculate and sum cooling rates (heat loss) due to interactions between
! electrons/neutrals, electrons/ions, ions/neutrals
!
      do i=lon0,lon1 ! (DO 500)
        do k=lev0,lev1-1
!
! Electron/N2 loss rate:
! loss_en2 = (L0(E,N2)+L0(E,N2,ROT)+L0(E,N2,VIB))/NE (s11)
!
        loss_en2(k,i) = n2n(k,i)*(1.77E-19*(1.-1.21E-4*te(k,i))*
     |    te(k,i) + 2.9e-14/root_te(k,i) + loss_en2v(k,i))
!
! Start total of electron/neutral loss rate (s11):
!
        loss_en(k,i) = loss_en2(k,i)
!
! Electron/O2 loss rates: (L0(E,O2)+L0(E,O2,ROT)+L0(E,O2,VIB)/NE
!
          loss_eo2(k,i) = o2n(k,i)*(1.21e-18*(1.+3.6e-2*root_te(k,i))*
     |      root_te(k,i)+6.9e-14/root_te(k,i)+3.125e-21*te(k,i)**2)
          loss_en(k,i) = loss_en(k,i)+loss_eo2(k,i)
!
! Electron/O(1d) loss rates: L0(E,O,1D)/(NE*N(O))
!
          loss_eo1d(k,i) = 22713.*(1./te(k,i)-1./tn(k,i))
          loss_eo1d(k,i) = sign(abs(loss_eo1d(k,i))+del,loss_eo1d(k,i))
          loss_eo1d(k,i) = 22713./(te(k,i)*tn(k,i))*
     |      (1.-exp(loss_eo1d(k,i)))/loss_eo1d(k,i)
!
! loss_eo1d function often fails here with bad argument to exp()
! due to high te and/or high loss_eo1d from above.
!         loss_eo1d(k,i) = 1.57e-12*exp((2.4e4+0.3*(te(k,i)-1500.)-
!    |      1.947e-5*(te(k,i)-1500.)*(te(k,i)-4000.))*(te(k,i)-3000.)/
!    |      (3000.*te(k,i)))*loss_eo1d(k,i)
          loss_eo1d(k,i) = 0.
!
! Electron/O1 loss rates: (L0(E,O)+L0(E,O,F))/NE
!
          loss_eo1(k,i) = o1n(k,i)*(7.9e-19*(1.+5.7e-4*te(k,i))*
     |      root_te(k,i)+3.4e-12*(1.-7.e-5*te(k,i))/tn(k,i)*
     |      (150./te(k,i)+0.4))

          loss_en(k,i) = loss_en(k,i)+loss_eo1(k,i)
!
! loss_xen = L0*(E,N) (s8)
!
          loss_xen(k,i) = (loss_en(k,i)+o1n(k,i)*(1.-alam/(ad+sd*
     |      n2n(k,i)))*loss_eo1d(k,i))*root_ne(k,i)*evergs
!
! Complete total electron/neutral loss rate L0(E,N) (s11):
!
          loss_en(k,i) = (loss_en(k,i)+o1n(k,i)*loss_eo1d(k,i))*
     |      root_ne(k,i)*evergs
!
! Calculate L0(E,I) = L(E,I)/(TE-TI), where L(E,I) is loss due to
! interactions between electrons and ions.
!
          loss_ei(k,i) = 3.2e-8*root_ne(k,i)/(root_te(k,i)*te(k,i))*
     |      15.*(op(k,i)+0.5*o2p(k,i)+0.53*nop(k,i))*evergs

          root_tn(k,i) = sqrt(tn(k,i))
! 
! loss_in = ion/neutral cooling = L0(I,N) =L(I,N)/(TI-TN)
!
          loss_in(k,i) = ((6.6e-14*n2n(k,i)+5.8e-14*o2n(k,i)+0.21e-14*
     |      o1n(k,i)*root2*root_tn(k,i))*op(k,i)+(5.45e-14*o2n(k,i)+
     |      5.9e-14*n2n(k,i)+4.5e-14*o1n(k,i))*nop(k,i)+(5.8e-14*
     |      n2n(k,i)+4.4e-14*o1n(k,i)+0.14e-14*o2n(k,i)*root_tn(k,i))*
     |      o2p(k,i))*evergs
!
! Complete tridiagonal matrix coefficients and rhs:
!
! q_coef = q_coef-(L0(E,N)+L0(E,I))/TE**2.5 = Q
!
          q_coef(k,i) = q_coef(k,i)-(loss_en(k,i)+loss_ei(k,i))/
     |      te(k,i)**2.5
!          
! rhs = rhs-L0(E,N)*TN-L0(E,I)*TI
!
          rhs(k,i) = rhs(k,i)-loss_en(k,i)*tn(k,i)-loss_ei(k,i)*ti(k,i)
! by Lei
          lossout(k,i) =  loss_en(k,i)*(te(k,i)-tn(k,i))
     |		        + loss_ei(k,i)*(te(k,i)-ti(k,i))

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1 (end DO 500)
!
! 8/28/06 btf: Jiuhou added EHEAT and ELOSS:
!     write(*,*)'come test tei',lossout(k,i), qout(k,i)
!     call addfsech('EHEAT'  ,' ',' ',qout    ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('ELOSS'  ,' ',' ',lossout ,lon0,lon1,nk,nkm1,lat)    

!     call addfsech('L_EN2'  ,' ',' ',loss_en2 ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_EO2'  ,' ',' ',loss_eo2 ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_EO1D' ,' ',' ',loss_eo1d,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_EO1'  ,' ',' ',loss_eo1 ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_XEN'  ,' ',' ',loss_xen ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_EN'   ,' ',' ',loss_en  ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_EI'   ,' ',' ',loss_ei  ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('L_IN'   ,' ',' ',loss_in  ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('Q_COEF' ,' ',' ',q_coef   ,lon0,lon1,nk,nkm1,lat)
!     call addfsech('RHS1'   ,' ',' ',rhs      ,lon0,lon1,nk,nkm1,lat)
!
! Calculate heating due to electron/neutral and electron/ion collisions
! (ergs/sec/gm):
!
      do i=lon0,lon1 ! (DO 24)
        do k=lev0,lev1-1
          if (te(k,i)-ti(k,i) >= 0.) then
            q_eni(k,i)=loss_ei(k,i)*(te(k,i)-ti(k,i))
          else
            q_eni(k,i) = 0.
          endif
          q_eni(k,i) = (loss_xen(k,i)*(te(k,i)-tn(k,i))+q_eni(k,i))
     |      *avo/xnmbar(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1 (end DO 24)

!     call addfsech('Q_ENI',' ',' ',q_eni,lon0,lon1,nk,nkm1,lat) ! s7
!
! Add collisional heating to Q for use in thermodynamic equation.
      do i=lon0,lon1 ! (DO 27)
        do k=lev0,lev1-2
          qtotal(k+1,i,lat) = qtotal(k+1,i,lat)+
     |      .5*(q_eni(k,i)+q_eni(k+1,i))
        enddo ! k=lev0,lev1-2
!
! Upper and lower boundaries:
        qtotal(lev0,i,lat) = qtotal(lev0,i,lat)+1.5*q_eni(lev0,i)-
     |    0.5*q_eni(lev0+1,i)
        qtotal(lev1,i,lat) = qtotal(lev1,i,lat)+1.5*q_eni(lev1-1,i)-
     |    0.5*q_eni(lev1-2,i)
      enddo ! i=lon0,lon1 ! (DO 27)

!     call addfsech('Q_TOT',' ',' ',qtotal(:,:,lat),lon0,lon1,
!    |  nk,nk,lat)
!
! Solve tridiagonal system:
!
!     call addfsech('P_COEF' ,' ',' ',p_coef,lon0,lon1,nk,nkm1,lat) ! s1
!     call addfsech('Q_COEF' ,' ',' ',q_coef,lon0,lon1,nk,nkm1,lat) ! s2
!     call addfsech('R_COEF' ,' ',' ',r_coef,lon0,lon1,nk,nkm1,lat) ! s3
!     call addfsech('RHS2'   ,' ',' ',rhs   ,lon0,lon1,nk,nkm1,lat) ! s4
!
!     subroutine trsolv(a,b,c,f,x,lev0,lev1,k1,k2,lon0,lon1,lonmax,lat,
!    |  idebug)

      call trsolv(p_coef,q_coef,r_coef,rhs,te_out(:,lon0:lon1),
     |  lev0,lev1,lev0,lev1-1,lon0,lon1,nlonp4,lat,1)
!
! Periodic points:
!     call periodic_f2d(te_out(:,lon0:lon1),lon0,lon1,lev1-lev0+1)

!     call addfsech('TE_SOLV',' ',' ',te_out(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!
! Te = Te**(2./7.):
      do i=lon0,lon1
        do k=lev0,lev1-1
          te_out(k,i) = te_out(k,i)**(2./7.)
        enddo
      enddo
!
! 10/21/03 btf: make this check after te*(2/7), rather than before.
!
! Te must be >= Tn:
      do i=lon0,lon1
        do k=lev0,lev1-1
          if (te_out(k,i) < tn(k,i)) te_out(k,i) = tn(k,i)
        enddo
      enddo
!
! 10/4/05: put spval in top level:
      te_out(lev1,:) = spval
!
!     call addfsech('TE_OUT',' ',' ',te_out(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!
! Set ion temperature output. Use joule heating qji_ti from sub 
! qjoule_ti (see qjoule.F). lev1 not calculated.
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          ti_out(k,i) = (qji_ti(k,i)*(xnmbar(k,i)/avo)+
     |      loss_ei(k,i)*te_out(k,i)+loss_in(k,i)*tn(k,i))/
     |      (loss_ei(k,i)+loss_in(k,i))
!
! ti must be at least as large as tn:
          if (ti_out(k,i) < tn(k,i)) ti_out(k,i) = tn(k,i)
        enddo
      enddo
!
! 10/4/05: put spval in top level:
      ti_out(lev1,:) = spval
!
!     call addfsech('TI_OUT',' ',' ',ti_out(:,lon0:lon1),
!    |  lon0,lon1,nk,nkm1,lat)
!
#ifdef VT
!     code = 126 ; state = 'settei' ; activity='ModelCode'
      call vtend(126,ier)
#endif
      end subroutine settei
