      program  teta_test
      implicit none
      integer, parameter   :: lm=50
      integer, parameter   :: lmp1=lm+1  
      real, parameter      :: g=9.80
      real, parameter      :: Rd=287.04
      real, parameter      :: Rv=461.5
      real, parameter      :: eps=Rd/Rv 
      real, parameter      :: cp=1004.6
      real, parameter      :: prf0=101325
      real, parameter      :: T0=288.
      real, parameter      :: gamma=0.0065
      real, parameter      :: rgog=rd*gamma/g 
      real, parameter      :: pt=2500.                      
      integer              :: i,j,l,k,ivi
      real                 :: pd,hgt
      real                 :: phis,phub,phlb,phub2,phlb2
      real                 :: dummy
      real,dimension(lm)   :: etal, deta,heta,teta,tveta1,teta1,tveta2,teta2,qeta
      real,dimension(lmp1) :: eta,peta,zeta
       
      do l=1,lm
         read(2000,*)k,eta(l),eta(l+1)
      enddo
      read(2002,*) pd,hgt
       
      do l=1,lmp1
         peta(l)=(prf0-pt)*eta(l)+pt 
	 zeta(l)=t0*(1.-((pt+eta(l)*(prf0-pt))/prf0)**rgog)/gamma
         print*,eta(l),peta(l),zeta(l)
      enddo

      do l=1,lm
         etal(l)=0.5*(eta(l)+eta(l+1))
      enddo
      do l=lm,1,-1
         read(2001,*) k,deta(l),etal(l)
      enddo
      do l=lm,1,-1
        read(2003,*) i,j,k,heta(l),teta(l),qeta(l)
      enddo


      phis=g*hgt
      phub=0.
      phub2=0.
      do ivi=1,lm
         l=lmp1-ivi
         phlb=phub
	 phlb2=phub2
         phub=g*heta(l)
	 phub2=zeta(l)
         tveta1(l)=-(phlb-phub)*(pt+etal(l)*pd)/(rd*deta(l)*pd)
         teta1(l)=tveta1(l)/(qeta(l)*0.61+1)

         tveta2(l)=(phub2-phlb2)*g/(rd*log(peta(l+1)/peta(l)))
!--------Virtual temperature, TV=T*qeta(l), Q is specific humidity;
!        (see pp. 63-65 in Fleagle & Businger, 1963)	 
	 teta2(l)=tveta2(l)/(qeta(l)*0.61+1)
	 teta2(l)=tveta2(l)/((1./eps-1)*qeta(l)+1)	 
	 write(6,'(i3,3f10.2,3f10.4)')l,zeta(l),peta(l),heta(l),teta(l),teta1(l),teta2(l)

      enddo

      end program teta_test
