#include "dims.h"
!
      SUBROUTINE ALTV
      use cons_module,only: len1,len2,kmax,kmaxp1,expz,p0,boltz,
     |  rmassinv_o2,rmassinv_o,rmassinv_n2
      implicit none
!
! This routine is called from dynamics (cb advnce) before rates.
! Output is TVIB(i,k) (in taskcommon /DSN2PAR/, file crates_tdep.h).
! (top level kmaxp1 of TVIB is *not* defined)
! TVIB is used in RK2 calculation in sub rates.
!
#include "params.h"
#include "fcom.h"
#include "index.h"
#include "buff.h"
#include "crates_tdep.h"
#include "phys.h"
!
! Local:
      real :: viba(5),vibb(5),vibc(5),vibd(5),vibe(5),
     |  alfa(zimxp,zkmxp),w0(zimxp),w00(zimxp), clam(zimxp,zkmxp),
     |  dl(zimxp,zkmxp),w1(zimxp,zkmxp),w2(zimxp,zkmxp),w(zimxp,zkmxp)
      real :: xno(zimxp,zkmxp),   ! number density of o
     +        expsi(zimxp,zkmxp)  ! exp(-z)
      integer :: i,k,ntk,npo2k,npo1k,kvib,ntek,nzk,nek,m,jj,j1
      real :: x1,cl,g,hn2,dn2,tdif,tvt,cl1,teta,a0
!
      viba=(/2.565,-6.525,-7.670,-9.174,-9.819/)
      vibb=(/8.782e-4,1.001e-2,1.092e-2,1.204e-2,1.243e-2/)
      vibc=(/2.954e-7,-3.066e-6,-3.369e-6,-3.732e-6,-3.850e-6/)
      vibd=(/-9.562e-11,4.436e-10,4.891e-10,5.431e-10,5.600e-10/)
      vibe=(/7.252e-15,-2.449e-14,-2.706e-14,-3.008e-14,-3.100e-14/)
!
! Convert O1 from mmr to cm3 (store in xno):
!
      do k=1,kmaxp1
        do i=1,len1
          expsi(i,k) = expz(k)
        enddo
      enddo
      ntk = nj+nt      ! tn
      npo2k = nj+nps   ! o2 (mmr)
      npo1k = nj+npo1  ! o  (mmr)
      do i=1,len2
        xno(i,1)  = f(i,npo1k)*
     +    (expsi(i,1)*p0/(boltz*f(i,ntk)*
     +    (f(i,npo2k)*rmassinv_o2+f(i,npo1k)*rmassinv_o+
     +    (1.-f(i,npo2k)-f(i,npo1k))*rmassinv_n2)))*rmassinv_o
      enddo
c
ccccccc  THE CALCULATION OF THE QUANTA FREQUENCY AND LAMDA  cccc
C       ****
C       ****   FOR TGCM22 USE
	KVIB=29
C       ****
C       ****
C       ****   FOR TGCM13 USE
C       KVIB=9
C       ****
C       ****
      ntek = nj+nte+kvib-2
      ntk = nj+nt+kvib-2
      nek = nj+ne+kvib-2
      nzk = nj+nz+kvib-2
      do k=kvib,kmax
        ntek = ntek+1
        ntk = ntk+1
        nek = nek+1
        nzk = nzk+1
        do i=1,len1 
	  w0(i)=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(i)=w0(i)+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(i)=(1.-exp(cl))*10.**x1
          endif
	  w(i,k)=w0(i)*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(i,k)=2.*sqrt(tdif/tvt)
	  dl(i,k)=sqrt(tdif*tvt)
        enddo ! i=1,len1
      enddo ! k=kvib,kmax
      do i=1,len1
        w0(i)=0.
        w00(i)=0.
      enddo
      do k=kvib+1,kmax
        jj=kmax+kvib-k
        j1=jj+1
	do i=1,len1
	x1=(clam(i,jj)-clam(i,j1))/2.
	cl=exp(-clam(i,jj))
	cl1=exp(-clam(i,j1))
	w0(i)=w0(i)+(w(i,jj)*cl+w(i,j1)*cl1)*x1
	w00(i)=w00(i)+(w(i,jj)/cl+w(i,j1)/cl1)*x1
	w1(i,jj)=w0(i)
	w2(i,jj)=w00(i)
	enddo
      enddo
      do i=1,len1
        w1(i,kmax)=0.
        w2(i,kmax)=0.
      enddo
      ntk = nj+nt+kvib-2
      do k=kvib,kmax
        ntk = ntk+1
        do i=1,len1
          x1=-3353./f(i,ntk)
          teta=exp(x1)
	  cl1=exp(clam(i,k))
          x1=1./cl1
	  a0=teta+dl(i,k)/clam(i,k)*(w0(i)*(cl1-x1)-cl1*w1(i,k)+x1
     1       *w2(i,k))
          if(a0.gt.0.3) a0=0.3
	  if(a0.lt.1.e-3) a0=1.e-3
!
! FINAL VIBRATIONAL QUANTA AND TEMPERATURE
!
	  alfa(i,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
      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