source: LMDZ6/branches/Amaury_dev/libf/phylmd/stratosphere_mask.F90 @ 5441

! $Id: stratosphere_mask.F90 5158 2024-08-02 12:12:03Z fhourdin $

SUBROUTINE stratosphere_mask(missing_val, pphis, t_seri, pplay, xlat)

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! determination of tropopause height and temperature from gridded temperature data

  ! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003)
  ! modified: 6/28/06 tjr
  ! adapted to LMDZ by C. Kleinschmitt (2016-02-15)
  ! committed to LMDz by O. Boucher (2016) with a mistake
  ! mistake corrected by O. Boucher (2017-12-11)

  ! input:  temp(nlon,nlat,nlev)  3D-temperature field
  !         ps(nlon,nlat)         2D-surface pressure field
  !         zs(nlon,nlat)         2D-surface height
  !         nlon                  grid points in x
  !         nlat                  grid points in y
  !         pfull(nlon,nlat,nlev) full pressure levels in Pa
  !         plimu                 upper limit for tropopause pressure
  !         pliml                 lower limit for tropopause pressure
  !         gamma                 tropopause criterion, e.g. -0.002 K/m

  ! output: p_tropopause(klon)    tropopause pressure in Pa with missing values
  !         t_tropopause(klon)    tropopause temperature in K with missing values
  !         z_tropopause(klon)    tropopause height in m with missing values
  !         stratomask            stratospheric mask withtout missing values
  !         ifil                  # of undetermined values

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE dimphy
  USE phys_local_var_mod, ONLY: stratomask
  USE phys_local_var_mod, ONLY: p_tropopause, z_tropopause, t_tropopause
  USE lmdz_print_control, ONLY: lunout, prt_level
  USE lmdz_yomcst

  IMPLICIT NONE

  REAL, INTENT(IN) :: missing_val ! missing value, also XIOS
  REAL, DIMENSION(klon), INTENT(IN) :: pphis   ! Geopotentiel de surface
  REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay   ! pression pour le mileu de chaque couche (en Pa)
  REAL, DIMENSION(klon), INTENT(IN) :: xlat    ! latitudes pour chaque point

  REAL, PARAMETER :: plimu = 45000.
  REAL, PARAMETER :: pliml = 7500.
  REAL, PARAMETER :: gamma = -0.002
  LOGICAL, PARAMETER :: dofill = .TRUE.
  REAL, DIMENSION(klon) :: tp
  REAL, DIMENSION(klev) :: t, p
  INTEGER :: i, k, ifil
  REAL :: ptrp, ttrp, ztrp, psrf, zsrf, pi

  pi = 4. * ATAN(1.)

  !--computing tropopause
  DO i = 1, klon
    DO k = 1, klev
      t(k) = t_seri(i, klev + 1 - k)
      p(k) = pplay(i, klev + 1 - k)
    ENDDO
    psrf = pplay(i, 1)
    zsrf = pphis(i) / RG           !--altitude de la surface
    CALL twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp)
    tp(i) = ptrp
    p_tropopause(i) = ptrp
    z_tropopause(i) = ztrp
    t_tropopause(i) = ttrp
  ENDDO

  !--filling holes in tp but not in p_tropopause
  IF (dofill) THEN
    ifil = 0
    DO i = 1, klon
      IF (ABS(tp(i) / missing_val - 1.0)<0.01) THEN
        !set missing values to very simple profile (neighbour averaging too expensive in LMDZ)
        tp(i) = 50000. - 20000. * cos(xlat(i) / 360. * 2. * pi)
        ifil = ifil + 1
      ENDIF
    ENDDO

  ENDIF

  DO i = 1, klon
    DO k = 1, klev
      IF (pplay(i, k)<tp(i)) THEN
        stratomask(i, k) = 1.0
      ELSE
        stratomask(i, k) = 0.0
      ENDIF
    ENDDO
  ENDDO

  IF (ifil>0 .AND. prt_level >5) THEN
    WRITE(lunout, *)'Tropopause: number of undetermined values =', ifil
  ENDIF

END SUBROUTINE stratosphere_mask

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! twmo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE twmo(missing_val, level, t, p, ps, zs, plimu, pliml, gamma, ptrp, ttrp, ztrp)

  ! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003)

  USE lmdz_yomcst

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: level
  REAL, INTENT(IN) :: missing_val
  REAL, INTENT(IN), DIMENSION(level) :: t, p
  REAL, INTENT(IN) :: plimu, pliml, gamma, ps, zs
  REAL, INTENT(OUT) :: ptrp, ttrp, ztrp

  REAL, parameter :: deltaz = 2000.0

  REAL :: faktor
  REAL :: pmk, pm, a, b, tm, dtdp, dtdz, dlnp, tdlnp
  REAL :: ag, bg, ptph, ttph, a0, b0
  REAL :: pm0, tm0, pmk0, dtdz0
  REAL :: p2km, asum, aquer
  REAL :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
  INTEGER :: icount, jj, j

  ptrp = missing_val
  ttrp = missing_val
  ztrp = missing_val

  faktor = -RG / RD

  DO j = level, 2, -1

    ! dt/dz
    pmk = 0.5 * (p(j - 1)**rkappa + p(j)**rkappa)        ! p**k at half level
    pm = pmk**(1. / rkappa)                    ! p at half level
    a = (t(j - 1) - t(j)) / (p(j - 1)**rkappa - p(j)**rkappa)    ! dT/dp^k
    b = t(j) - (a * p(j)**rkappa)
    tm = a * pmk + b                    ! T at half level
    dtdp = a * rkappa * pm**(rkappa - 1.) ! dT/dp at half level
    dtdz = faktor * dtdp * pm / tm            ! dT/dz at half level

    ! dt/dz valid?
    IF (j==level)     go to 999                 ! no, start level, initialize first
    IF (dtdz<=gamma)  go to 999                 ! no, dt/dz < -2 K/km
    IF (dtdz0>gamma) go to 999                 ! no, dt/dz below > -2 K/km
    IF (pm>plimu)    go to 999                 ! no, pm too low

    ! dtdz is valid, calculate tropopause pressure ptph
    ag = (dtdz - dtdz0) / (pmk - pmk0)
    bg = dtdz0 - (ag * pmk0)
    ptph = exp(log((gamma - bg) / ag) / rkappa)

    ! calculate temperature at this ptph assuming linear gamma
    ! interpolation
    ttph = tm0
    ttph = ttph - (bg * log(pm0) + ag * (pm0**rkappa) / rkappa) / faktor * t(j)
    ttph = ttph + (bg * log(ptph) + ag * (ptph**rkappa) / rkappa) / faktor * t(j)

    IF (ptph<pliml) go to 999
    IF (ptph>plimu) go to 999

    ! 2nd test: dtdz above 2 km must not exceed gamma
    p2km = ptph + deltaz * (pm / tm) * faktor      ! p at ptph + 2km
    asum = 0.0                               ! dtdz above
    icount = 0                               ! number of levels above

    ! test until apm < p2km
    DO jj = j, 2, -1

      pmk2 = .5 * (p(jj - 1)**rkappa + p(jj)**rkappa)    ! p mean ^kappa
      pm2 = pmk2**(1 / rkappa)                      ! p mean
      IF(pm2>ptph) go to 110                ! doesn't happen
      IF(pm2<p2km) go to 888                ! ptropo is valid

      a2 = (t(jj - 1) - t(jj))                     ! a
      a2 = a2 / (p(jj - 1)**rkappa - p(jj)**rkappa)
      b2 = t(jj) - (a2 * p(jj)**rkappa)               ! b
      tm2 = a2 * pmk2 + b2                     ! T mean
      dtdp2 = a2 * rkappa * (pm2**(rkappa - 1))        ! dt/dp
      dtdz2 = faktor * dtdp2 * pm2 / tm2
      asum = asum + dtdz2
      icount = icount + 1
      aquer = asum / float(icount)               ! dt/dz mean

      ! discard ptropo ?
      IF (aquer<=gamma) go to 999           ! dt/dz above < gamma

      110 continue
    enddo                   ! test next level

    888 continue                    ! ptph is valid
    ptrp = ptph
    ttrp = ttph

    ! now calculate height of tropopause by integrating hypsometric equation
    ! from ps to ptrp
    ! linearly interpolate in p (results are identical to p^kappa)

    jj = LEVEL                                       ! bottom
    DO while ((P(jj)>PS) .OR. (T(jj)<100))     ! T must be valid too
      jj = jj - 1
    END DO

    DLNP = log(PS / P(jj))                             ! from surface pressure
    TM = T(jj)                                       ! take TM of lowest level (better: extrapolate)
    TDLNP = TM * DLNP

    DO while ((JJ>=2) .AND. (PTRP<P(jj - 1)))
      DLNP = log(P(jj) / P(jj - 1))
      TM = 0.5 * (T(jj) + T(jj - 1))
      TDLNP = TDLNP + TM * DLNP
      JJ = JJ - 1
    END DO

    DLNP = log(P(jj) / PTRP)                           ! up to tropopause pressure
    TM = 0.5 * (T(jj) + TTRP)                        ! use TTRP to get TM of this level
    TDLNP = TDLNP + TM * DLNP

    ZTRP = ZS + TDLNP * RD / RG

    !!if (ZTRP .lt. 0) THEN
    !!  PRINT*,'ZTRP=',ZTRP
    !!  PRINT*,'PS=',PS
    !!  PRINT*,'P=',P
    !!  PRINT*,'T=',T
    !!  PRINT*,'ZS=',ZS
    !!  stop
    !!endif

    return

    999 continue                    ! continue search at next higher level
    tm0 = tm
    pm0 = pm
    pmk0 = pmk
    dtdz0 = dtdz
    a0 = a
    b0 = b

  END DO

  ! no tropopouse found

END SUBROUTINE  twmo