!
      SUBROUTINE ALTV
!
! This routine is called from advnce before rates.
! Output is TVIB(i,k) (in common /RATES/, file crates.h).
! (top level kmaxp1 of TVIB is *not* defined)
! TVIB is used in RK2 calculation in sub rates.
!
      include 'params.h'
      include 'blnk.h'
      include 'index.h'
      include 'buff.h'
      include 'cons.h'
      include 'crates.h' ! this routine defines tvib(i,k)
!     include 'phys.h'   ! need j only for addfsech
!
      dimension viba(5),vibb(5),vibc(5),vibd(5),vibe(5)
      dimension alfa(zkmx)
      dimension clam(200),dl(200),w1(200),w2(200),w(200)
      data viba/2.565,-6.525,-7.670,-9.174,-9.819/
      data vibb/8.782e-4,1.001e-2,1.092e-2,1.204e-2,1.243e-2/
      data vibc/2.954e-7,-3.066e-6,-3.369e-6,-3.732e-6,-3.850e-6/
      data vibd/-9.562e-11,4.436e-10,4.891e-10,5.431e-10,5.600e-10/
      data vibe/7.252e-15,-2.449e-14,-2.706e-14,-3.008e-14,-3.100e-14/
      real :: xno(zimxp,zkmxp),   ! number density of o
     +        expsi(zimxp,zkmxp)  ! exp(-z)
!
! Convert O1 from mmr to cm3 (store in xno):
!
      do k=1,kmax
        expsi(:,k) = exps(k)
      enddo
      ntk = nj+nt      ! tn
      npo2k = nj+nps   ! o2 (mmr)
      npo1k = nj+nps2  ! o  (mmr)
      do i=1,len2
        xno(i,1)  = f(i,npo1k)*
     +    (expsi(i,1)*c(81) / (c(84)*f(i,ntk)*
     +    (f(i,npo2k)/rmass(1)+f(i,npo1k)/rmass(2)+
     +    (1.-f(i,npo2k)-f(i,npo1k))/rmass(3))))/rmass(2)
      enddo
c
ccccccc  THE CALCULATION OF THE QUANTA FREQUENCY AND LAMDA  cccc
C       ****
C       ****   FOR TGCM22 USE
C       KVIB=29
C       ****
C       ****
C       ****   FOR TGCM13 USE
      KVIB=9
C       ****
C       ****
      ntek = nj+nte+kvib-1
      ntk = nj+nt+kvib-1
      nek = nj+ne+kvib-1
      nzk = nj+nz+kvib-1
      do k=kvib,kmax
        ntek = ntek+1
        ntk = ntk+1
        nek = nek+1
        nzk = nzk+1
        do i=1,len1 
          w0=0.
          if(f(i,ntek) >= 1500.) then
            do m=1,5
              x1=viba(m)+f(i,ntek)*(vibb(m)+f(i,ntek)*(vibc(m)+vibd(m)*
     +          f(i,ntek)+vibe(m)*f(i,ntek)*f(i,ntek)))-16.
              cl=(3353./f(i,ntek)-3353./f(i,ntk))*float(m)
              x1=(1.-exp(cl))*10.**x1
              w0=w0+x1
            enddo
          else
            x1=-5.922+f(i,ntek)*(3.151e-2+f(i,ntek)*(-4.075e-5+
     +        f(i,ntek)*2.439e-8-f(i,ntek)*f(i,ntek)*5.479e-12))-16.
            cl=3353./f(i,ntek)-3353./f(i,ntk)
            w0=(1.-exp(cl))*10.**x1
          endif
          w(k)=w0*f(i,nek)
          g=981./(1.+f(i,nzk)*1.e-5/6378.)**2 ! 1.e-5 for Z in km
          hn2=826.e5*f(i,ntk)/g/28.
          x1=xno(i,k)/1.e9
          dn2=9.69e7*f(i,ntk)**0.724/x1
          tdif=hn2*hn2/dn2
          tvt=1./(0.107*exp(-69.9/f(i,ntk)**0.33)*x1)
          clam(k)=2.*sqrt(tdif/tvt)
          dl(k)=sqrt(tdif*tvt)
        enddo ! i=1,len1
      enddo ! k=kvib,kmax
!
      w0=0.
      w00=0.
      do k=kvib+1,kmax
        jj=kmax+kvib-k
        j1=jj+1
        x1=(clam(jj)-clam(j1))/2.
        cl=exp(-clam(jj))
        cl1=exp(-clam(j1))
        w0=w0+(w(jj)*cl+w(j1)*cl1)*x1
        w00=w00+(w(jj)/cl+w(j1)/cl1)*x1
        w1(jj)=w0
        w2(jj)=w00
      enddo
      w1(kmax)=0.
      w2(kmax)=0.
!
      ntk = nj+nt+kvib-1
      do k=kvib,kmax
        ntk = ntk+1
        do i=1,len1
          x1=-3353./f(i,ntk)
          teta=exp(x1)
          cl1=exp(clam(k))
          x1=1./cl1
          a0=teta+dl(k)/clam(k)*(w0*(cl1-x1)-cl1*w1(k)+x1*w2(k))
          if(a0.gt.0.3) a0=0.3
!
! FINAL VIBRATIONAL QUANTA AND TEMPERATURE
!
          alfa(k)=a0
          tvib(i,k)=-3353./log(a0/(1.+a0))
!         if (tvib(i,k).lt.f(i,ntk)) tvib(i,k)=f(i,ntk)
        enddo
      enddo
      ntk = nj+nt-1
      do k=1,kvib-1
        ntk = ntk+1
        do i=1,len1
          tvib(i,k)=f(i,ntk)
        enddo
      enddo
!     call addfsech('TVIB',tvib,zimxp,zkmxp,kmax,j)
      return
      end
c
cccccccccccccccccccccc SEPTEMBER 1997 VERSION cccccccccccccccc
c     N2 VIBRATIONAL QUANTA STEADY STATE
c     APPROXIMATION is given by equation (B.3) :
c
c     Pavlov, A.V., The role of vibrationally excited nitrogen in the
c     formation of the mid-latitude negative ionospheric storms,
c     Annales Geophysicae, 1994, Vol.12, No. 10, P. 554-564.
c
c     Accuracy of the steady state approximation is evaluated by
c
c     Pavlov, A.V. and M.J.Buonsanto, Using steady state vibrational
c     temperatures to model effects of N2* on calculations of electron
c     densities, J. Geophys. Res., Vol.101, No. A12, P.26941-26945, 1996.
c
c     This subroutine uses the revised constants for calculations of the
c     production frequency of N2 vibrational quanta by thermal electron
c     excitation of N2 given by
c
c     Pavlov, A.V., New electron energy transfer rates for vibrational
c     excitation of N2, Annales Geophysicae, 1997, Vol.15,
c     in press (accepted).
c
c     INPUT:
c     h(nn)  - ALTITUDE(KM)
c
c     h(1)<140 km !!!  h(nn)>400 km !!!!  h(n)-h(n-1)<20 km !!!
c
c     yo(nn) - O NUMBER DENSITY(CM-3)
c     tn(nn) - NEUTRAL TEMPERATURE (K)
c     te(nn) - ELECTRON TEMPERATURE (K)
c     cne(nn)- ELECTRON NUMBER DENSITY (CM-3)
c     nn - a number of altitude points
c
c    !!!!!  nn<200    !!!!!
c
c     OUTPUT:
c     alfa(nn) - N2 VIBRATIONAL QUANTA
c     tv(nn)   - N2 VIBRATIONAL TEMPERATURE  (K)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc