pro find_trop,z,p,t,zt,pt,tt

;---------------------------------------------------------------
;  Find the tropopause altitude from the input profiles
;
; Input:
;  z.......altitude (km), vector
;  p.......pressure, vector
;  t.......temperature (K), vector
;
; Output:
;  zt......tropopause altitude (km)
;  pt......tropopause presure
;  tt......tropopause temperature (K)
;
; Source: Mark Hervig
;---------------------------------------------------------------

      nl   = n_elements(z)
      pt   = 0.0
      zt   = 0.0
      tt   = 0.0
      zmax = 20.
      zmin = 7.
      nlay = 7  ; find dt/dz in this many layers

;---- call the trop where dT/dz becomes less than dtdz_trop (K/km),
;     i.e. look for a notable slope change,  set to work bottom up

      if (z(1) gt z(2)) then begin
        n1 = nl-1
        n2 = nlay
        n3 = -1
        n4 = -1*nlay
      endif 
      if (z(1) lt z(2)) then begin
        n1 = 0
        n2 = nl-nlay-2 
        n3 = 1
        n4 = nlay
      endif 
 
      for i = n1,n2,n3 do begin
        if (z(i) le zmax and z(i) ge zmin) then begin

          sxx  = 0.0
          syy  = 0.0
          sxy  = 0.0
          sxsq = 0.0
          sysq = 0.0

          for j = i,i+n4,n3 do begin
            xax  = z(j)-z(i)
            yax  = t(j)
            sxsq = sxsq+xax^2
            sysq = sysq+yax^2
            sxy  = sxy+xax*yax
            sxx  = sxx+xax
            syy  = syy+yax
          endfor
          sxy = sxy-syy*sxx/float(nlay+1)
          sxx = sxsq-sxx*sxx/float(nlay+1)
          syy = sysq-syy*syy/float(nlay+1)

          if (sxx eq 0.0) then begin
            beta = 1.e25
          endif else begin
            beta = sxy/sxx
          endelse

          if (beta gt -2.0) then begin
            tt = t(i)
            pt = p(i)
            zt = z(i)
            goto,jump1
          endif

        endif
      endfor

      jump1:
      return
      end