MODULE tropopause_m

!  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
  PRIVATE
  PUBLIC :: dyn_tropopause

CONTAINS

!-------------------------------------------------------------------------------
!
SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
!
!-------------------------------------------------------------------------------
  USE assert_m,      ONLY: assert
  USE assert_eq_m,   ONLY: assert_eq
  USE comvert_mod,   ONLY: preff
  USE dimphy,        ONLY: klon, klev
  USE geometry_mod,  ONLY: latitude_deg, longitude_deg
  USE interpolation, ONLY: locate
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
  REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
  REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
  REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
  REAL, INTENT(IN)  :: presnivs(:) !--- Cells-centers nominal pressure
  REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
  REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
  REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
!-------------------------------------------------------------------------------
! Local variables:
  include "YOMCST.h"
  CHARACTER(LEN=80)  :: sub
  INTEGER :: it, i, nsrf, k, nh, kmin, kmax
  REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
  REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
  REAL, DIMENSION(klon,klev-1) :: pvort
!-------------------------------------------------------------------------------
  sub='dyn_tropopause'
  CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
  CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
  CALL assert(SHAPE(rot)  ==[klon,klev  ],TRIM(sub)//" rot shape")

  !--- DEFAULT THRESHOLDS
  tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
  pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
  pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
  pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
  kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
  kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1

  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
  DO k= 1, klev
    DO i = 1, klon
      t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
    END DO
  END DO

  !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
  t_edg(:,1) = ts(:)
  DO k= 2, klev
    DO i = 1, klon
      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
      t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
    END DO
  END DO

  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
  DO k= 1, klev-1
    DO i = 1, klon
      IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
      pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
                         * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
                           - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
                         / ( paprs(i,k+1)  -     paprs(i,k  ) )
    END DO
  END DO

  !--- LOCATE TROPOPAUSE
  ptrop(:)=0.
  DO i = 1, klon
    !--- Dynamical tropopause
    it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
    IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
      .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
      al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
    END IF
    !--- Potential temperature iso-surface
    it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
    IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
      al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
    END IF
  END DO

END SUBROUTINE dyn_tropopause

END MODULE tropopause_m