var/da/da_physics/da_trop_wmo.inc

   1 subroutine da_trop_wmo (t, z, p, nlev, tropt, tropp, tropk)
   2
   3    !----------------------------------------------------------------------------
   4    ! * Computes tropopause T, P, and/or level based on code from Cameron Homeyer
   5    !     and WMO definition
   6    !
   7    ! * WMO tropopause definition:
   8    !     The boundary between the troposphere and the stratosphere, where an
   9    !     abrupt change in lapse rate usually occurs. It is defined as the lowest
  10    !     level at which the lapse rate decreases to 2 °C/km or less, provided
  11    !     that the average lapse rate between this level and all higher levels
  12    !     within 2 km does not exceed 2 °C/km.
  13    !----------------------------------------------------------------------------
  14
  15    implicit none
  16
  17    ! Assumed shape inputs for single column (size=nlev)
  18    !  ordered from bottom to top of model
  19    real, intent(in) :: t(:)             ! Temperature, K. (3D)
  20
  21 !JJG: "z" is supposed to be height, not geopotential height. Does it matter?
  22    real, intent(in) :: z(:)             ! Geopotential height above m.s.l., m.
  23    real, intent(in) :: p(:)             ! Pressure, mb.
  24    real, optional, intent(out) :: tropt ! Tropopause temperature, K.
  25    real, optional, intent(out) :: tropp ! Tropopause pressure, mb.
  26    integer, optional, intent(out) :: tropk
  27    integer, intent(in) :: nlev
  28
  29    real    :: dtdz,laps  !LAPS
  30    integer :: dtdztest(nlev), ztest(nlev)
  31    integer :: i, j, k, kk, ktrop
  32
  33 !   real, parameter :: tropz_min = 5000.0
  34 !   real, parameter :: tropz_max = 19000.0
  35
  36    if (.not.present(tropt) .and. &
  37        .not.present(tropp) .and. &
  38        .not.present(tropk)) return
  39
  40    if (present(tropt)) tropt = missing_r
  41    if (present(tropp)) tropp = missing_r
  42
  43    !Loop over levels to find tropopause (single column)
  44    ktrop = nlev-1
  45    trop_loop: do  k = 1, nlev-1
  46       if ( p(k) .le. 500.0 ) then
  47          ! Compute lapse rate (-dT/dz)
  48          dtdz = ( t(k+1) - t(k)   ) / &
  49                 ( z(k)   - z(k+1) )
  50       else
  51          ! Set lapse rate for p > 500 hPa
  52          dtdz = 999.9
  53       endif
  54       !Check if local lapse rate <= 2 K/km
  55       if (dtdz .le. 0.002) then
  56          ! Initialize lapse rate and altitude test arrays
  57          dtdztest = 0
  58          ztest = 0
  59
  60          ! Compute average lapse rate across levels above current candidate
  61          do kk = k+1, nlev-1
  62             dtdz = ( t(kk+1) - t(k)   ) / &
  63                    ( z(k)    - z(kk+1) )
  64
  65             !If avg. lapse rate <= 2 K/km and z <= trop + 2 km, set pass flag
  66             if ( ( dtdz .le. 0.002 ) .and. &
  67                  ( (z(k) - z(kk)) .le. 2000. ) ) THEN
  68                dtdztest(kk) = 1
  69             endif
  70
  71             ! If z <= trop + 2 km, set pass flag
  72             IF ( (z(k) - z(kk)) .le. 2000.0 ) THEN
  73                ztest(kk) = 1
  74             endif
  75          enddo   !kk loop
  76          laps=dtdz !LAPS
  77          IF (SUM(dtdztest) .eq. SUM(ztest)) THEN
  78             ! If qualified as tropopause, set altitude index and return value
  79             ktrop = k
  80             exit trop_loop
  81          ENDIF
  82       ENDIF
  83    end do trop_loop
  84
  85 !   ! Filter ktrop using tpause height thresholds
  86 !   ztest = 0
  87 !   if ( z(ktrop) .gt. tropz_max ) then
  88 !      where ( z.le.tropz_max )
  89 !         ztest = 1
  90 !      end where
  91 !      do  k = nlev, 1, -1
  92 !         if (ztest(k) .eq. 1) then
  93 !            ktrop = k
  94 !            exit
  95 !         end if
  96 !      end do
  97 !   else if ( z(ktrop) .lt. tropz_min ) then
  98 !      where ( z.ge.tropz_min )
  99 !         ztest = 1
 100 !      end where
 101 !      do  k = 1, nlev
 102 !         if (ztest(k) .eq. 1) then
 103 !            ktrop = k
 104 !            exit
 105 !         end if
 106 !      end do
 107 !   end if
 108
 109    if (present(tropt)) tropt = t(ktrop)
 110    if (present(tropp)) tropp = p(ktrop)
 111    if (present(tropk)) tropk = ktrop
 112
 113 end subroutine da_trop_wmo
 114