! Version 1.0 released by David Romps on September 12, 2017.
! Version 1.1 vectorized lcl.R, released on May 24, 2021.
! 
! When using this code, please cite:
! 
! @article{16lcl,
!   Title   = {Exact expression for the lifting condensation level},
!   Author  = {David M. Romps},
!   Journal = {Journal of the Atmospheric Sciences},
!   Year    = {2017},
!   Month   = dec,
!   Number  = {12},
!   Pages   = {3891--3900},
!   Volume  = {74}
! }
!
! This lcl function returns the height of the lifting condensation level
! (LCL) in meters.  The inputs are:
! - p in Pascals
! - T in Kelvins
! - Exactly one of rh, rhl, and rhs (dimensionless, from 0 to 1):
!    * The value of rh is interpreted to be the relative humidity with
!      respect to liquid water if T >= 273.15 K and with respect to ice if
!      T < 273.15 K. 
!    * The value of rhl is interpreted to be the relative humidity with
!      respect to liquid water
!    * The value of rhs is interpreted to be the relative humidity with
!      respect to ice
! - return_ldl is an optional logical flag.  If true, the lifting deposition
!   level (LDL) is returned instead of the LCL. 
! - return_min_lcl_ldl is an optional logical flag.  If true, the minimum of the
!   LCL and LDL is returned.

program test

   use lcl_mod

   implicit none

   double precision :: p, rh, T1, T2

   p = 1.e5
   rh = 0.5
   T1 = 300.
   T2 = 200.

   if ( &
      abs(lcl(p,T1,rhl=rh,return_ldl=.false.)/( 1433.844139279)-1.) < 1.e-10 .and. &
      abs(lcl(p,T1,rhs=rh,return_ldl=.false.)/( 923.2222457185)-1.) < 1.e-10 .and. &
      abs(lcl(p,T2,rhl=rh,return_ldl=.false.)/( 542.8017712435)-1.) < 1.e-10 .and. &
      abs(lcl(p,T2,rhs=rh,return_ldl=.false.)/( 1061.585301941)-1.) < 1.e-10 .and. &
      abs(lcl(p,T1,rhl=rh,return_ldl=.true. )/( 1639.249726127)-1.) < 1.e-10 .and. &
      abs(lcl(p,T1,rhs=rh,return_ldl=.true. )/( 1217.336637217)-1.) < 1.e-10 .and. &
      abs(lcl(p,T2,rhl=rh,return_ldl=.true. )/(-8.609834216556)-1.) < 1.e-10 .and. &
      abs(lcl(p,T2,rhs=rh,return_ldl=.true. )/( 508.6366558898)-1.) < 1.e-10 ) then
      print *,'Success'
   else
      print *,'Failure'
   end if

end program test