***********************************************************************
*     This program reads a radiosonde profile and interpolates it to a
*     finer resolution vertical array, so that there will typically be
*     no more than 1 K or so temperature change between interpolated
*     levels.
*
*     The input consists of series of levels, starting at the surface.
*     Each level is given on a separate line and consists of pressure
*     (mb), temperature (C), and dewpoint (C).
*
*     The output has the same format.
*
*     Typically, this program is compiled and executed as follows (under
*     Unix/Linux):
*
*     f77 -o interp_sounding interp_sounding.f 
*     interp_sounding < inputfile > outputfile
*
*     G. Petty (March 10, 2002)
************************************************************************

      parameter (NMAX = 100, NNMAX = 1000)
      real p(NMAX), t(NMAX), td(NMAX)
      real pp(NNMAX), tt(NNMAX), ttd(NNMAX)

      n = 1
 1    read(*,*,err=888,end=999)p(n),t(n),td(n)

* Check for valid order
      if (n .ge. 2) then
         if (p(n-1) .le. p(n)) then
            stop 'Levels must be entered starting with highest pressure'
         endif
      endif
* Check for valid values
      if (p(n) .le. 0.0) stop 'Invalid negative or zero pressure'
      if (p(n) .gt. 1050.0) stop 'Pressure must be given in millibars'
      if (t(n) .gt. 40.0) stop 'Temperature must be in Celsius'
      if (td(n) .gt. t(n)) stop 'Invalid dewpoint temperature'

      if (n .lt. NMAX) then
         n = n + 1
         goto 1
      else
         stop 'Too many levels in input'
      endif
 888  continue
      stop 'Format error in input'
      
 999  continue

*     actual number of levels successfully read
      nlevels = n - 1

*     If we're here, we've successfully read in a valid sounding.  Now start
*     the postprocessing.

*     Interpolate sounding to approximately 0.15 km intervals

*     Variable ZZ is non-dimensional pseudo-altitude, measured in
*     pressure scale heights.

      p0 = p(1)
      zz1 = 0.0
      nn = 1
* copy surface level in sounding to interpolated profile
      pp(1) = p(1)
      tt(1) = t(1)
      ttd(1) = td(1)

* interpolate between pairs of levels in original sounding
      do n=1,nlevels-1
         zz2 = -log(p(n+1)/p0)
         delzz = zz2 - zz1
         intervals = nint(delzz/0.02)
         if (intervals .eq. 0) intervals = 1
         dzz = delzz/intervals

         delt = t(n+1)-t(n)
         deltd = td(n+1)-td(n)

         if (delt .lt. -5.0 .and. delt/delzz .lt. -100.0) then
            pause
     &           'Superadiabatic lapse rate detected - check input!'
         endif
         do i=1,intervals
            nn = nn + 1
            if (nn .gt. NNMAX) 
     &           stop 'Too many interpolated levels - fix NNMAX'
            
            zznew = i*dzz 
            pp(nn) = p0*exp(-(zz1 + zznew))
            tt(nn) = t(n) + zznew*delt/delzz
            ttd(nn) = td(n) + zznew*deltd/delzz

         enddo
         zz1 = zz2
      enddo
      nnlevels = nn

* output interpolated profile

      do nn=1,nnlevels
         write(*,100)pp(nn),tt(nn),ttd(nn)
 100     format(f10.5,f8.3,f8.3)
      enddo

      stop
      end