subroutine stattrop(pfull, tfull, nlev, ptrop, itrop)

c**********************************************************************************
c
c     input : p, t        real(nlev)
c     input : nlev        real
c     output: ptrop       real
c     output: itrop     integer
c
c     programmed by Dominik Brunner              V1.0 Aug 2000
c
c     built upon routine stattrop by Peter van Velthoven,
c     KNMI, The Netherlands
c     and on the ECHAM4/CHEM routine xttphwmo by
c     Thomas Reichle, Christine Land, B. Steil and R. Hein, DLR
c
c     purpose
c     -------
c     stattrop computes the pressure (Pa) at the thermal (static) tropopause
c     (TP) for a 1D column (sounding) of pressures and temperatures following
c     the definition of the height of the tropopause as postulated by WMO (1957).
c
c     ATTENTION: In the current formulation of the program the first
c     level (index 1) must be at the top of the atmosphere and the
c     last level (index nlev) at the surface
c
c     interface
c     ---------
c     call stattrop(pfull, tfull, nlev, ptrop, itrop)
c     - Input
c       nlev : number of vertical levels
c       pfull: pressure in 1D column at nlev full levels (layers)
c       tfull: temperature in 1D column at nlev full levels
c     - Output
c       ptrop: height of the tropopause in Pa
c       itrop: index of layer containing the tropopause
c
c     methods
c     -------
c     - Lapse rate gamma = -dT/dz
c       Using the hydrostatic approximation
c
c       dz = -dp/(g*rho) = -dp/p * R/g * T = -dlnp * R/g * T
c
c       we get -dT/dz = dT/T * g/R *1/dlnp = dlnT/dlnp
c
c     - Variables are assumed to vary linearly in log(pressure)
c     - The tropopause is the lowest level above 450 hPa fullfilling
c       the WMO criterium. If ptrop is < 85 hPa it is set to 85 hPa.
c       If no tropopause is found ptrop is set to -999.
c
c     references
c     ----------
c
c     - WMO (1992): International meteorological vocabulary, Genf, 784pp.:
c
c       " 1. The first tropopause is defined as the lowest level at which
c       the lapse rate decreases to 2 deg C per kilometer or less,
c       provided also the average lapse rate between this level and
c       all higher levels within 2 kilometers does not exceed 2 deg C"
c     
c     - Randel WJ, Wu F, Gaffen DJ, Interannual variability of the tropical 
c       tropopause derived from radiosonde data and NCEP reanalyses,
c       JOURNAL OF GEOPHYSICAL RESEARCH, 105, 15509-15523, 2000.
c
c       The following webpage clearifies the calculation of the tropopause
c       in the NCEP reanalysis: http://dss.ucar.edu/pub/reanalysis/FAQ.html
c
c       For comparison NCEP reanalysis tropopause pressures can be obtained
c       from http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml
c
c**********************************************************************************

      implicit none
      integer nlev
      real pfull(nlev), tfull(nlev)
      real ptrop, p2km, lapseavg, lapsesum
      integer ilev, ilev1, itrop, count, l
      integer nlevm
      parameter (nlevm=31)
      real lapse(nlevm), phalf(nlevm)
      real grav, rgas, const, lapsec, rinval
      parameter (grav=9.80665, rgas=287.04, const=1000.*grav/rgas, 
     &     lapsec=2.0, rinval=-999.0 )
c     min. and max. allowed tropopause pressure
      real pmin, pmax
      parameter (pmin=8500.,pmax=45000.)
c     deltaz is layer thickness used for second tropopause criterium
      real deltaz
      parameter (deltaz=2.)
      logical found
      real p1, p2 , weight
      
      ptrop = rinval
      itrop = 0
      found = .false.

c     Loop from model top to surface and calculate lapse rate gamma=-dT/dz (K/km)
      do ilev = 1, nlev-1
         ilev1 = ilev + 1
         lapse(ilev) = ALOG(tfull(ilev))- ALOG(tfull(ilev1))
         lapse(ilev) = lapse(ilev) /
     &        (ALOG(pfull(ilev))- ALOG(pfull(ilev1)))
         lapse(ilev) = const * lapse(ilev)
         phalf(ilev) = (pfull(ilev) + pfull(ilev1))*0.5
      end do
c     Loop from surface to top to find (lowest) tropopause
      do ilev = nlev-2, 1, -1
c     **** 1st test: -dT/dz less than 2 K/km and pressure LT pmax? ****
         if (lapse(ilev).lt.lapsec.and.pfull(ilev).lt.pmax) then
            if (.not.found) then
c     Interpolate tropopause pressure linearly in log(p)
               p1     = ALOG(phalf(ilev))
               p2     = ALOG(phalf(ilev+1))
               if (lapse(ilev).ne.lapse(ilev+1)) then
                  weight = (lapsec - lapse(ilev)) /
     &                 (lapse(ilev+1)-lapse(ilev))
                  ptrop = exp(p1 + weight * (p2-p1)) ! tropopause pressure
               else
                  ptrop = phalf(ilev)
               end if
c     *** 2nd test: average -dT/dz in layer deltaz above TP must not be greater than 2K/km
               p2km = exp(ALOG(ptrop)-deltaz*const/tfull(ilev)) ! press 2 km above TP
               lapsesum = 0.0  ! init avg. lapse rate 2 km above TP
               lapseavg=1.e9
               count = 0 
               do l=ilev,1,-1
                  if (phalf(l).lt.p2km) goto 100
                  lapsesum=lapsesum+lapse(l)
                  count=count+1
                  lapseavg=lapsesum/float(count)
               end do
 100           continue
               found=lapseavg.lt.lapsec
c     Discard tropopause?
               if (.not.found) then
                  ptrop=rinval
               else
                  ptrop=max(ptrop,pmin) ! ptrop must be GE 85 hPa
                  itrop=ilev            ! assign level index
               end if
            end if
         end if
      end do
c      if (.not.found) then
c         write(*,*)
c     &        'STATTROP-warning: no static tropopause found:'
c         write(*,*) 'phalf:', phalf
c         write(*,*) 'lapse:', lapse
c      end if
c      
      return
      end