subroutine thetae(tk, p, rh, pl, tl, w, te)
! This routine calculates the equivalent potential temperature
! using the formula described in Bolton (1980) MWR
!
! Input: tk  - Temperature (K)
!        p   - Total pressure (hPa)
!        rh  - Relative humidity (%)
!
! Output: tl - Temperature (K) of the LCL
!         pl - Pressure (hPa) of the LCL
!         te - Thetae (K)
!         w  - Mixing ratio (g/Kg)

! Note: If any input values is outside the normal
!       atmospheric range, all output variables are set to -999.

! Modifications:
! 4/20/21 (J. Dostalek) Started this Fortran 90 version
!-------------------------------------------------------------------------------

  implicit none

  real, intent(in) :: tk
  real, intent(in) :: p
  real, intent(inout) :: rh

  real, intent(out) :: tl
  real, intent(out) :: pl
  real, intent(out) :: te
  real, intent(out) :: w

  real, parameter :: ctk = 273.15

  real :: t
  real :: es
  real :: ws
  real :: aa
  real :: bb

  logical :: endsub

! Check input
  endsub = .false.
  if (tk < 100. .or. tk >  350.) endsub = .true.
  if (p  <=  0. .or. p  > 1200.) endsub = .true.
  if (rh <  -5. .or. rh >  105.) endsub = .true.

  if (endsub) then
    pl = -999.
    tl = -999.
    te = -999.
    w  = -999.
  else
    if (rh <=   0.) rh =  1.0
    if (rh >= 100.) rh = 99.9

    t = tk - ctk

    es = 6.112*exp(17.67*t/(t + 243.5))
    ws = 622.0*es/(p - es)
    w = ws*rh/100.0

    tl = 55.0 + 1.0/(1.0/(tk - 55.0) - alog(rh/100.0)/2840.)
    pl = p*(tl/tk)**(1.0/0.2854)

    aa = .2854*(1.0 - .00028*w)
    bb = ((3.376/tl) - .00254)*w*(1.+.00081*w)
    te = (tk*(1000.0/p)**aa)*exp(bb)
  end if

  return
end subroutine thetae