!
      subroutine settei(tn,o2,ox,n2,ne,op,o2p,nop,barm,qji_ti,
     |  te,ti,te_out,ti_out,lev0,lev1,lon0,lon1,lat)
!
! Calculate electron and ion temperatures.
!
      use params_module,only: dz,nlonp4
      use cons_module,only: pi,rtd,evergs,expzmid_inv,expzmid,
     |  expz,rmassinv_o2,rmassinv_o1,rmassinv_n2,boltz,p0,grav,
     |  gask,dipmin,avo
      use input_module,only: f107,step,step_compqrj
      use qrj_module,only:
     |  qtotal, ! total heating             ! F(NQ)
     |  qo2p,   ! o2+ ionization            ! F(NQO2P)
     |  qop,    ! o+  ionization            ! F(NQOP)
     |  qn2p,   ! n2+ ionization            ! F(NQN2P)
     |  qnop,   ! no+ ionization            ! F(NQNOP)
     |  qnp,    ! n+  ionization            ! F(NQNP)
     |  qop2p,  ! o+(2p) ionization         ! F(NQOP2P)
     |  qop2d   ! o+(2d) ionization         ! F(NQOP2D)
      use magfield_module,only: rlatm,dipmag
      use chapman_module,only: chi   ! solar zenith angle (nlonp4,nlat)
      use init_module,only: istep ! istep for debug only
      use addfld_module,only: addfld
      implicit none
!
! Input args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,       ! neutral temperature    (f(:,nj+nt  ))
     |  o2,       ! molecular oxygen       (f(:,nj+nps ))
     |  ox,       ! ox (mmr)               (f(:,nj+nps2))
     |  n2,       ! n2 (mmr)               (f(:,nj+npn2))
     |  ne,       ! electron density (cm3) (f(:,nj+ne  ))
     |  op,       ! o+                     (f(:,nj+nop ))
     |  o2p,      ! o2+                    (f(:,nj+no2p))
     |  nop,      ! no+                    (f(:,nj+nnop))
     |  barm,     ! mean molecular weight  (f(:,nj+nms ))
     |  qji_ti    ! joule heating from qjoule_ti (qjoule.F)
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  te_out,  ! output electron temperature (deg K) (f(:,njnp+nte))
     |  ti_out   ! output ion temperature (deg K)      (f(:,njnp_nti))
!
! Inout args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) ::
     |  te,ti    ! from previous step
!
! Local:
      integer :: i,k,n,i0,i1
      real,parameter ::
     |  fpolar = -5.0e+8,  ! 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)
!
! Quantities at interfaces:
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  tei,       ! te at interfaces (full levels) ! s14
     |  o2i,       ! o2 at interfaces (full levels) ! s13
     |  oxi,       ! ox at interfaces (full levels) ! s12
     |  tni,       ! tn at interfaces (full levels) ! s10
     |  n2i,       ! n2 at interfaces (full levels) ! s11
     |  o2n,       ! O2 number density (interfaces)
     |  oxn,       ! OX number density (interfaces)
     |  n2n,       ! N2 number density (interfaces)
     |  xnmbari,   ! p0*e(-z)*barm/kT  (interfaces)
     |  xnmbarm,   ! p0*e(-z)*barm/kT  (midpoints)
     |  root_te,   ! sqrt(te)
     |  root_tn,   ! sqrt(tn)
     |  root_ne,   ! sqrt(ne)                    (s15)
     |  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)
!
! 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)
!
! For diagnostic plotting (addfld):
      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
     |  sindipmag_ki ! for diagnostic plotting of sindipmag
!
      real,dimension(lev0:lev1) :: expzp ! exp(-z)
!
      i0=lon0 ; i1=lon1
      root2 = sqrt(2.)
      f107te = f107
      if (f107te > 235.) f107te = 235.

!     call addfld('QJI_TI',' ',' ',qji_ti(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
      do i=lon0,lon1
        if (abs(rlatm(i,lat))-pi/4.5 >= 0.) then
          a(i) = 1.                       ! t3
        else
          a(i) = .5*(1.+sin(pi*(abs(rlatm(i,lat))-pi/9.)/(pi/4.5)))
        endif 
!
! Increased heat flux for TE fom protonosphere.
!!      fed(i) = ( -5.0e+7*f107te*a(i)-4.0e+7*f107te)*1.2 ! t4
	fed(i) = ( -5.0e+7*f107te*a(i)-4.0e+7*f107te)     ! t4
        fen(i) = fed(i)/2.                ! t5
!       fed(i) = fed(i)+qteaur(i,lat)     ! tiegcm only
!       fen(i) = fen(i)+qteaur(i,lat)     ! tiegcm only
        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)
        fed_ki(:,i) = fed(i)
        fen_ki(:,i) = fen(i)
        fe_ki (:,i) = fe(i)
      enddo ! i=lon0,lon1

!     call addfld('MAGLAT',' ',' ',a_ki     ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('CHI'   ,' ',' ',chi_ki   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FED'   ,' ',' ',fed_ki   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FEN'   ,' ',' ',fen_ki   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('FE'    ,' ',' ',fe_ki    ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

!
! Set limits on te,ti from previous step (timegcm only, not in tiegcm):
! Write warning to stdout if limits are exceeded:
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          if (te(k,i) < tn(k,i)) then
!           write(6,"('>>> settei: lat=',i3,' k=',i3,' i=',i3,
!    |        ' te < tn: te=',e12.4,' setting te=tn=',e12.4)") 
!    |        lat,k,i,te(k,i),tn(k,i)
            te(k,i) = tn(k,i)
          endif
          if (ti(k,i) < tn(k,i)) then
!           write(6,"('>>> settei: lat=',i3,' k=',i3,' i=',i3,
!    |        ' ti < tn: te=',e12.4,' setting ti=tn=',e12.4)") 
!    |        lat,k,i,ti(k,i),tn(k,i)
            ti(k,i) = tn(k,i)
          endif
          if (te(k,i) > 5000.) then
!           write(6,"('>>> settei: lat=',i3,' k=',i3,' i=',i3,
!    |        ' te > 5000: te=',e12.4,' setting te=5000.')") 
!    |        lat,k,i,te(k,i)
            te(k,i) = 5000.
          endif
          if (ti(k,i) > 5000.) then
!           write(6,"('>>> settei: lat=',i3,' k=',i3,' i=',i3,
!    |        ' ti > 5000: ti=',e12.4,' setting ti=5000.')") 
!    |        lat,k,i,ti(k,i)
            ti(k,i) = 5000.
          endif
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
! 
! te,o2,o,tn at interfaces: 
!
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          tei(k,i) = .5*(te(k,i)+te(k-1,i)) ! s14
          o2i(k,i) = .5*(o2(k,i)+o2(k-1,i)) ! s13
          oxi(k,i) = .5*(ox(k,i)+ox(k-1,i)) ! s12
          tni(k,i) = .5*(tn(k,i)+tn(k-1,i)) ! s10
        enddo ! k=lev0+1,lev1-1
!
! Bottom:
        tei(lev0,i) = 1.5*te(lev0,i)-.5*te(lev0+1,i)
        o2i(lev0,i) = 1.5*o2(lev0,i)-.5*o2(lev0+1,i)
        oxi(lev0,i) = 1.5*ox(lev0,i)-.5*ox(lev0+1,i)
        tni(lev0,i) = 1.5*tn(lev0,i)-.5*tn(lev0+1,i)
!
! Top:
        tei(lev1,i) = 1.5*te(lev1-1,i)-.5*te(lev1-2,i)
        o2i(lev1,i) = 1.5*o2(lev1-1,i)-.5*o2(lev1-2,i)
        oxi(lev1,i) = 1.5*ox(lev1-1,i)-.5*ox(lev1-2,i)
        tni(lev1,i) = 1.5*tn(lev1-1,i)-.5*tn(lev1-2,i)
!
! N2:
        do k=lev0,lev1
          n2i(k,i) = 1.-oxi(k,i)-o2i(k,i)        ! s11
          if (n2i(k,i) < 0.) n2i(k,i) = 1.e-20   ! s11
          if (o2i(k,i) < 0.) o2i(k,i) = 1.e-20   ! s13
          if (tei(k,i) < 100.) tei(k,i) = 100.   ! s14
          if (tni(k,i) < 100.) tni(k,i) = 100.   ! s10
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfld('TEI',' ',' ',tei,'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O2I',' ',' ',o2i,'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OXI',' ',' ',oxi,'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('TNI',' ',' ',tni,'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('N2I',' ',' ',n2i,'ilev',lev0,lev1,'lon',i0,i1,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
          xnmbari(k,i) = p0*expzp(k)*barm(k,i)/(boltz*tni(k,i)) ! s8
          o2i(k,i) = xnmbari(k,i)*o2i(k,i)*rmassinv_o2       ! s13
          oxi(k,i) = xnmbari(k,i)*oxi(k,i)*rmassinv_o1       ! s12
          n2i(k,i) = xnmbari(k,i)*n2i(k,i)*rmassinv_n2       ! s11
          root_te(k,i) = sqrt(tei(k,i))
!
          tek0(k,i) = 7.5e5/
     |      (1.+3.22e4*tei(k,i)**2/ne(k,i)*(root_te(k,i)*
     |      (2.82e-17 - 3.41e-21   * tei    (k,i))*n2i(k,i)+
     |      (2.20e-16 + 7.92e-18   * root_te(k,i))*o2i(k,i)+
     |       1.10e-16 * (1.+5.7e-4 * tei (k,i))*oxi(k,i)))*evergs
!
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfld('XNMBARI',' ',' ',xnmbari,
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('TEK0'   ,' ',' ',tek0   ,
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!
      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*tni(k,i)/(barm(k,i)*grav)                 ! s6
        enddo ! k=lev0,lev1
      enddo ! i=lon0,lon1

!     call addfld('H_MID'   ,' ',' ',h_mid(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('H_INT'   ,' ',' ',h_int,
!    |  'ilev',lev0,lev1,'lon',i0,i1,lat)

      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 ! tiegcm only
!
! 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)*tni(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 addfld('DIPMAG' ,' ',' ',dipmag_ki,
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('SINDIPM',' ',' ',sindipmag_ki,
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('P_COEF' ,' ',' ',p_coef(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('Q_COEF' ,' ',' ',q_coef(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('R_COEF' ,' ',' ',r_coef(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('RHS0'   ,' ',' ',rhs(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

!
! qtot = total ionization rate = sum(Qxx) = 
! (QO2+) + (QO+) + (QN2+) + (QNO+) + (QN+)
!
      call addfld('QO2P' ,' ',' ',qo2p (:,:,lat),
     |  'ilev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('QOP'  ,' ',' ',qop  (:,:,lat),
     |  'ilev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('QN2P' ,' ',' ',qn2p (:,:,lat),
     |  'ilev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('QNOP' ,' ',' ',qnop (:,:,lat),
     |  'ilev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('QNP'  ,' ',' ',qnp  (:,:,lat),
     |  'ilev',lev0,lev1,'lon',i0,i1,lat)
!
! Note re comparison with tgcm24: Because the check for qtot < 1.e-20
! is inside the reduction operation loop, timegcm1 must sum separately
! for each ion species, as in tgcm24. If timegcm1 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,5
        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)
        end select
        do i=lon0,lon1
          do k=lev0,lev1
            qtot(k,i) = qtot(k,i)+fki(k,i)                       ! s11
            if (qtot(k,i) < 1.e-20) qtot(k,i) = 1.e-20
          enddo
        enddo
      enddo ! n=1,7

!     call addfld('QTOT_SUM',' ',' ',qtot,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)

      do i=lon0,lon1
        do k=lev0,lev1-1
          qtot(k,i) = sqrt(qtot(k,i)*qtot(k+1,i))+              ! s11 (timegcm only)
     |                sqrt(qop2p(k,i,lat)*qop2p(k+1,i,lat))+
     |                sqrt(qop2d(k,i,lat)*qop2d(k+1,i,lat))
          root_ne(k,i) = sqrt(ne(k,i)*ne(k+1,i))                ! s15
          xnmbarm(k,i) = p0*expz(k)*.5*(barm(k,i)+barm(k+1,i))/ ! s6
     |      (boltz*tn(k,i))
          o2n(k,i) = xnmbarm(k,i)*o2(k,i)*rmassinv_o2           ! s14
          oxn(k,i) = xnmbarm(k,i)*ox(k,i)*rmassinv_o1           ! s13
          n2n(k,i) = xnmbarm(k,i)*n2(k,i)*rmassinv_n2
!
! 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*oxn(k,i))) ! s10
          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.3
!
! Subtract qe from right-hand-side:
          rhs(k,i) = rhs(k,i)-qe(k,i)*qtot(k,i)*evergs          ! s4
          root_te(k,i) = sqrt(te(k,i))                          ! s5
!
! 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))         ! s9
          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))        ! s10
          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
        qtot(lev1,i) = 0.
      enddo ! i=lon0,lon1

!     call addfld('QTOT_SUM',' ',' ',qtot   ,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('ROOT_NE' ,' ',' ',root_ne(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('XNMBARM' ,' ',' ',xnmbarm(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('O2N'     ,' ',' ',o2n(lev0:lev1-1,:)    ,
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('OXN'     ,' ',' ',oxn(lev0:lev1-1,:)    ,
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('N2N'     ,' ',' ',n2n(lev0:lev1-1,:)    ,
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('QE'      ,' ',' ',qe(lev0:lev1-1,:)     ,
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('RHS'     ,' ',' ',rhs(lev0:lev1-1,:)    ,
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('ROOT_TE' ,' ',' ',root_te,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('C_EN2V'  ,' ',' ',coll_en2v(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('L_EN2V'  ,' ',' ',loss_en2v(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

!
! Calculate and sum cooling rates (heat loss) due to interactions between
! electrons/neutrals, electrons/ions, ions/neutrals
!
      do i=lon0,lon1
        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))*    ! s11
     |      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)                                ! s11
!
! 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))             ! s10
          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.
! Gpi and/or ncep runs may stop here when TE gets too hot, e.g., >
!   8000-10000 K or too small (< 10.) It often gets too hot around
!   lat=-60, slt=20.
! (loss_eo1d is zeroed out in some versions of tiegcm1)
!
! 5/21/04 btf: setting loss_eo1d==0. for timegcm1:
! 5/25/04 btf: using the loss_eo1d function again:
!
          loss_eo1d(k,i) = 1.57e-12*exp((2.4e4+0.3*(te(k,i)-1500.)-  ! s10
     |      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.                                        ! tiegcm1
!
! Electron/O1 loss rates: (L0(E,O)+L0(E,O,F))/NE
!
          loss_eo1(k,i) = oxn(k,i)*(7.9e-19*(1.+5.7e-4*te(k,i))*     ! s11
     |      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)+oxn(k,i)*(1.-alam/(ad+sd*    ! s8
     |      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)+oxn(k,i)*loss_eo1d(k,i))*     ! s11
     |      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))* ! s10
     |      15.*(op(k,i)+0.5*o2p(k,i)+0.53*nop(k,i))*evergs

          root_tn(k,i) = sqrt(tn(k,i))                               ! s5
! 
! 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* ! s9
     |      oxn(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*oxn(k,i))*nop(k,i)+(5.8e-14*
     |      n2n(k,i)+4.4e-14*oxn(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))/     ! s2
     |      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) ! s4

        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

      call addfld('L_EN2'  ,' ',' ',loss_en2(lev0:lev1-1,:) ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_EO2'  ,' ',' ',loss_eo2(lev0:lev1-1,:) ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_EO1D' ,' ',' ',loss_eo1d(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_EO1'  ,' ',' ',loss_eo1(lev0:lev1-1,:) ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_XEN'  ,' ',' ',loss_xen(lev0:lev1-1,:) ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_EN'   ,' ',' ',loss_en(lev0:lev1-1,:)  ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_EI'   ,' ',' ',loss_ei(lev0:lev1-1,:)  ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('L_IN'   ,' ',' ',loss_in(lev0:lev1-1,:)  ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('Q_COEF' ,' ',' ',q_coef(lev0:lev1-1,:)   ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('RHS1'   ,' ',' ',rhs(lev0:lev1-1,:)      ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,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)) 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/xnmbarm(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1 (end DO 24)

!     call addfld('Q_ENI' ,' ',' ',q_eni(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! Add collisional heating to Q for use in thermodynamic equation.
!
! 1/2/08 btf: Increment qtotal only if qrj has been called.
!   (note qtotal is not used in this routine, only incremented
!    by q_eni for later use in dt)
!
      if (istep==1.or.mod(istep,step_compqrj/step)==0) then
        do i=lon0,lon1
          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
      endif

!     call addfld('QTOT_TEI',' ',' ',qtotal(:,:,lat),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Solve tridiagonal system:
!
      call addfld('P_COEF' ,' ',' ',p_coef(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('Q_COEF' ,' ',' ',q_coef(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('R_COEF' ,' ',' ',r_coef(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('RHS2'   ,' ',' ',rhs(lev0:lev1-1,:)   ,
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
!     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 addfld('TE_SOLV',' ',' ',te_out(:,lon0:lon1),
     |  'lev',lev0,lev1,'lon',lon0,lon1,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
!
! Te must be >= Tn:
      do i=lon0,lon1
        do k=lev0,lev1-1
          if (te_out(k,i) < tn(k,i)) then
            te_out(k,i) = tn(k,i)
          endif
        enddo
      enddo
!
! Te is not defined at lev1 (only up to lev1-1)
      call addfld('TE_OUT',' ',' ',te_out(:,lon0:lon1),
     |  'lev',lev0,lev1,'lon',lon0,lon1,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)*(xnmbarm(k,i)/avo)+
     |      loss_ei(k,i)*te(k,i)+loss_in(k,i)*tn(k,i))/  ! tiegcm uses te_out(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
!     call addfld('QJI_TI',' ',' ',qji_ti(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!     call addfld('TI_OUT',' ',' ',ti_out(:,lon0:lon1),
!    |  'lev',lev0,lev1,'lon',lon0,lon1,lat)
!
#ifdef VT
!     code = 126 ; state = 'settei' ; activity='ModelCode'
      call vtend(126,ier)
#endif
      end subroutine settei
