subroutine tfromzp()
c
c calculate the temperature profile based upon the
c height as a function of pressure.
c
c input
c  via commons:
c plginv - log of pressure levels (11)
c z11 - altitude of pressure levels
c nx - # of levels
c r0 - radius of earth (km)
c tmptop=tmpprf(1) - temperature of top layer (k)
c p101lg - log of 101 pressure levels
c
c  on call line
c tmpsfc - temperature of surface (k)
c 
c output:
c t_101 - avg. temperature for the 101 layers
c n_101 - # of layers for t_101
c
c
      implicit none
      include 'parameter.inc'
      include 'atm_11.cmn'
      include 'input.cmn'
      include 'atm_101.cmn'
      include 'consts.cmn'
      include 'switches.cmn'
c local
      integer mode,i
      real*8 z1,z2,const,t2,gamma,tmptop
c
      mode=0
      tmptop=tmpprf(1)
      const=1.0/(-34.163)
      if(jprint(6).ne.0)then
        write(33,'(a9)')'tfromzp. '
        write(33,'(a40)')'p_101lg(i) exp(p_101lg(i)) z2 t_101(i)'
      endif
c
c compute the average temperature of each of the 101 pressure layers
c
c spline has already been set
          call splint2(plginv,z11,cp,n_p_lay,plginv(n_p_lay),z2,t2,mode)
          t2=t2/const/((1.0+z2/r)**2)
          gamma=2.*(tmptop-t2)/(81.-z11(n_p_lay))*tmptop*const*
     &((1.0+z11(n_p_lay)/r)**2)
          z1=0.0
          do i=n_101,1,-1
             call splint2(plginv,z11,cp,n_p_lay,p_101lg(i),z2,t_101(i),
     &mode)
             if(p_101lg(i).lt.plginv(n_p_lay))then
c use constant lapse rate (as a function of ln(pres)
c determined by the temperature at plginv(n_p_lay)
c and the average value of the temperature between z11(n_p_lay) and 81 km
c which should be tmptop.
               t_101(i)=t2+gamma*(p_101lg(i)-plginv(n_p_lay))
cc*****evaluate weighted temperatures in the upper most 40 layers
cc     based on constant lapse rate of -1.5 degrees per layer)
c               t_101(i)=t_101(i+1)-1.5
               z2=z1+0.5*(t_101(i+1)+t_101(i))*const*((1.0+z1/r)**2)*
     &(p_101lg(i)-p_101lg(i+1))
             else
c use the slope of z(ln(p)) and correct for gravity variation
               t_101(i)=t_101(i)/const/((1.0+z2/r)**2)
             endif
             alt_101(i)=z2
             if(jprint(6).ne.0)
     &      write(33,'(1P4e12.4)')p_101lg(i),exp(p_101lg(i)),
     &z2,t_101(i)
             z1=z2
          enddo
c now average the values
          do i=101,2,-1
             t_101(i)=0.5*(t_101(i)+t_101(i-1))
          enddo
      return
      end