!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_ES                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_es( p, t ) 7,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation specific humidity using enhanced Teten's
!  formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/08/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_es     Saturation water vapor pressure (Pa)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_es      ! Saturation water vapor pressure (Pa)
!
!-----------------------------------------------------------------------
!
!  Function f_es and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_esl, f_esi

!fpp$ expand (f_esl)
!fpp$ expand (f_esi)
!dir$ inline always f_esl, f_esi
!*$*  inline routine (f_esl, f_esi)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( t >= 273.15 ) THEN      ! for water
    f_es = f_esl( p,t )
  ELSE                            ! for ice
    f_es = f_esi( p,t )
  END IF

  RETURN
  END FUNCTION f_es

!
!-----------------------------------------------------------------------
!
!  Calculate the saturation water vapor over liquid water using
!  enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_esl( p, t ) 3

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_esl     ! Saturation water vapor pressure over liquid water

  REAL :: f

  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  f = satfwa + satfwb * p
  f_esl = f * satewa * EXP( satewb*(t-273.15)/(t-satewc) )

  RETURN
  END FUNCTION f_esl
!
!-----------------------------------------------------------------------
!
!  Calculate the saturation water vapor over ice using enhanced
!  Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_esi( p, t ) 3

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_esi     ! Saturation water vapor pressure over ice (Pa)

  REAL :: f

  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  f = satfia + satfib * p
  f_esi = f * sateia * EXP( sateib*(t-273.15)/(t-sateic) )

  RETURN
  END FUNCTION f_esi
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_DESDT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_desdt( t ) 2,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate d(es)/dt/es using enhanced Teten's formula. See function
!  f_es.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/12/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_desdt  d(es)/dt/es (1/K)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: t         ! Temperature (K)
  REAL :: f_desdt   ! d(es)/dt/es (1/K)
!
!-----------------------------------------------------------------------
!
!  Function f_desdtl and f_desdti and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_desdtl, f_desdti

!fpp$ expand (f_desdtl)
!fpp$ expand (f_desdti)
!dir$ inline always f_desdtl, f_desdti
!*$*  inline routine (f_desdtl, f_desdti)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( t >= 273.15 ) THEN      ! for water
    f_desdt = f_desdtl(t)
  ELSE                         ! for ice
    f_desdt = f_desdti(t)
  END IF

  RETURN
  END FUNCTION f_desdt
!
!-----------------------------------------------------------------------
!
!  Calculate d(esl)/dt/esl using enhanced Teten's formula.
!  See function f_esl.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_desdtl( t ) 7

  IMPLICIT NONE

  REAL :: t         ! Temperature (K)
  REAL :: f_desdtl  ! d(esl)/dt/esl (1/K)

  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  f_desdtl = satewb*(273.15-satewc)/(t-satewc)**2

  RETURN
  END FUNCTION f_desdtl
!
!-----------------------------------------------------------------------
!
!  Calculate d(esi)/dt/esi using enhanced Teten's formula.
!  See function f_esi.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_desdti( t ) 4

  IMPLICIT NONE

  REAL :: t         ! Temperature (K)
  REAL :: f_desdti  ! d(esi)/dt/esi (1/K)

  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  f_desdti = sateib*(273.15-sateic)/(t-sateic)**2

  RETURN
  END FUNCTION f_desdti
!
!-----------------------------------------------------------------------
!
!  Calculate d(qvs)/dt/qvs using enhanced Teten's formula.
!  See function f_es.
!
!           Rd
!  rddrv = ----
!           Rv
!
!  qvs = rddrv * es / (p - ( 1.0 - rddrv ) * es )
!
!    1    d(qvs)           1 - rddrv             1    d(es)
!  ----- -------- = ( 1 + ----------- * qvs ) * ---- -------
!   qvs     dT               rddrv               es    dT
!
!-----------------------------------------------------------------------
!


  FUNCTION f_dqvsdt( p, t ),2

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pa)
  REAL :: t         ! Temperature (K)
  REAL :: f_dqvsdt  ! d(esi)/dt/esi (1/K)

  REAL :: qvsat, desdt
!
!-----------------------------------------------------------------------
!
!  Function f_desdtl and f_desdti and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_desdt, f_qvsat

!fpp$ expand (f_desdt)
!fpp$ expand (f_qvsat)
!dir$ inline always f_desdt, f_qvsat
!*$*  inline routine (f_desdt, f_qvsat)

  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  qvsat = f_qvsat( p,t )
  desdt = f_desdt( t )
  f_dqvsdt = ( 1. + qvsat*(1.-rddrv)/rddrv ) * desdt

  RETURN
  END FUNCTION f_dqvsdt
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_QVSAT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_qvsat( p, t ) 38,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation specific humidity using enhanced Teten's
!  formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/08/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_qvsat  Saturation water vapor specific humidity (kg/kg).
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_qvsat   ! Saturation water vapor specific humidity (kg/kg)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_es and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_es

!fpp$ expand (f_es)
!dir$ inline always f_es
!*$*  inline routine (f_es)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  f_qvsat = rddrv * f_es(p,t) / (p-(1.0-rddrv)*f_es(p,t))

  RETURN
  END FUNCTION f_qvsat
!
!-----------------------------------------------------------------------
!
!  Calculate the saturation specific humidity over liquid water using
!  enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_qvsatl( p, t ) 13,1

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_qvsatl  ! Saturation specific humidity over liquid
                    ! water (kg/kg)

  REAL :: fesl

  INCLUDE 'phycst.inc'

  REAL :: f_esl

!fpp$ expand (f_esl)
!dir$ inline always f_esl
!*$*  inline routine (f_esl)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  fesl = f_esl( p, t )
  f_qvsatl = rddrv * fesl / (p-(1.0-rddrv)*fesl)

  RETURN
  END FUNCTION f_qvsatl
!
!-----------------------------------------------------------------------
!
!  Calculate the saturation specific humidity over ice using
!  enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_qvsati( p, t ) 5,1

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_qvsati  ! Saturation specific humidity over ice (kg/kg)

  REAL :: fesi

  INCLUDE 'phycst.inc'

  REAL :: f_esi

!fpp$ expand (f_esi)
!dir$ inline always f_esi
!*$*  inline routine (f_esi)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  fesi = f_esi( p, t )
  f_qvsati = rddrv * fesi / (p-(1.0-rddrv)*fesi)

  RETURN
  END FUNCTION f_qvsati
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_MRSAT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_mrsat( p, t ) 9,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation water vapor mixing ratio using enhanced
!  Teten's formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/08/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    f_mrsat  Saturation water vapor mixing ratio (kg/kg).
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_mrsat   ! Saturation water vapor mixing ratio (kg/kg)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: fes
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_es and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_es

!fpp$ expand (f_es)
!dir$ inline always f_es
!*$*  inline routine (f_es)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  fes = f_es( p,t )
  f_mrsat = rddrv * fes / (p-fes)

  RETURN
  END FUNCTION f_mrsat
!
!-----------------------------------------------------------------------
!
!  Calculate the saturated water vapor mixing ratio over liquid water
!  using enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_mrsatl( p, t ),1

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_mrsatl  ! Saturation water vapor mixing ratio over liquid
                    ! water (kg/kg)

  REAL :: fesl

  INCLUDE 'phycst.inc'

  REAL :: f_esl

!fpp$ expand (f_esl)
!dir$ inline always f_esl
!*$*  inline routine (f_esl)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  fesl = f_esl( p,t )
  f_mrsatl = rddrv * fesl / (p-fesl)

  RETURN
  END FUNCTION f_mrsatl
!
!-----------------------------------------------------------------------
!
!  Calculate the saturated water vapor mixing ratio over ice
!  using enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_mrsati( p, t ),1

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: f_mrsati  ! Saturation water vapor mixing ratio over ice
                    ! (kg/kg)

  REAL :: fesi

  INCLUDE 'phycst.inc'

  REAL :: f_esi

!fpp$ expand (f_esi)
!dir$ inline always f_esi
!*$*  inline routine (f_esi)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  fesi = f_esi( p,t )
  f_mrsati = rddrv * fesi / (p-fesi)

  RETURN
  END FUNCTION f_mrsati
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_TDEW                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_tdew( p, t, qv ) 2,3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the dew point temperature using reversed enhanced
!  Tetan's formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/09/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!    qv       Specific humidity (kg/kg)
!
!  OUTPUT:
!
!    f_tdew   Dew point temperature (K)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: t         ! Temperature (K)
  REAL :: qv        ! Specific humidity (kg/kg)
  REAL :: f_tdew    ! Dew point temperature (K)
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_tdewl and f_tdewi and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_tdewl, f_tdewi

!fpp$ expand (f_tdewl)
!fpp$ expand (f_tdewi)
!dir$ inline always f_tdewl, f_tdewi
!*$*  inline routine (f_tdewl, f_tdewi)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( t >= 273.15 ) THEN         ! for water
    f_tdew = f_tdewl( p,qv )
    IF ( f_tdew < 273.15 ) THEN  ! using ice formula
      f_tdew = f_tdewi( p,qv )
    END IF
  ELSE                            ! for ice
    f_tdew = f_tdewi( p,qv )
  END IF

  RETURN
  END FUNCTION f_tdew
!
!-----------------------------------------------------------------------
!
!  Calculate the dew point temperature over liquid water using the
!  reversed enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_tdewl( p, qv ) 1

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: qv        ! Specific humidity (kg/kg)
  REAL :: f_tdewl   ! Dew point temperature over liquid water (K)

  REAL :: qvs, esl, f, ln

  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  qvs = MAX( qv, 1.e-8 )

  f = satfwa + satfwb * p
  esl = p*qvs/(rddrv+(1-rddrv)*qvs)
  ln = ALOG( esl/(f*satewa) )

  f_tdewl = ( ln*satewc - 273.15*satewb ) / ( ln - satewb )

  RETURN
  END FUNCTION f_tdewl
!
!-----------------------------------------------------------------------
!
!  Calculate the dew point temperature over ice using the
!  reversed enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!


  FUNCTION f_tdewi( p, qv ) 2

  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascal)
  REAL :: qv        ! Specific humidity (kg/kg)
  REAL :: f_tdewi   ! Dew point temperature over ice (K)

  REAL :: qvs, esi, f, ln

  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  qvs = MAX( qv, 1.e-8 )

  f = satfia + satfib * p
  esi = p*qvs/(rddrv+(1-rddrv)*qvs)
  ln = ALOG( esi/(f*sateia) )

  f_tdewi = ( ln*sateic - 273.15*sateib ) / ( ln - sateib )

  RETURN
  END FUNCTION f_tdewi
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE F_PT2PTE                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_pt2pte( p, pt, qv ) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the equivalent potential temperature of an air
!  parcel that is either saturated or unsaturated.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/14/1998  Re-created from subroutine EQUIPT
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    p        Base state pressure (Pascals)
!    pt       Potential temperature (degrees Kelvin)
!    qv       Water vapor specific humidity (kg/kg)
!
!  OUTPUT:
!
!    f_pt2pte Equivalent potential temperature (K)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p         ! Pressure (Pascals)
  REAL :: pt        ! Potential temperature (degrees Kelvin)
  REAL :: qv        ! Water vapor specific humidity (kg/kg)

  REAL :: f_pt2pte  ! Equivalent potential temperature (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: t       ! Temperature
  REAL :: tdew    ! Dew point temperature
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_tdew and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_tdew

!fpp$ expand (f_tdew)
!dir$ inline always f_tdew
!*$*  inline routine (f_tdew)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  t = pt * (p/p0)**rddcp

  tdew = f_tdew( p,t,qv )    ! Dew point temperture.

  f_pt2pte = pt * EXP( lathv*qv/(cp*tdew) )

  RETURN
  END FUNCTION f_pt2pte
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETQVS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getqvs( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend,             & 14,1
           p, t, qvsat )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation specific humidity using Teten's formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/17/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/03/92 (M. Xue)
!  Further documentation.
!
!  2/10/93 (K. Droegemeier)
!  Cleaned up documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!  12/11/1997 (Yuhe Liu)
!  Rewrote the subroutine calling function F_QVSAT to calculate the
!  saturated specific humidity for array.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Array dimension in x-direction
!    ny       Array dimension in y-direction
!    nz       Array dimension in z-direction
!    ibgn     Starting index in x-direction
!    jbgn     Starting index in y-direction
!    kbgn     Starting index in z-direction
!    iend     Last index in x-direction
!    jend     Last index in y-direction
!    kend     Last index in z-direction
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    qvsat    Saturation water vapor specific humidity (kg/kg).
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! Input array dimension

  INTEGER :: ibgn             ! Starting index in x-direction
  INTEGER :: jbgn             ! Starting index in y-direction
  INTEGER :: kbgn             ! Starting index in z-direction
  INTEGER :: iend             ! Last index in x-direction
  INTEGER :: jend             ! Last index in y-direction
  INTEGER :: kend             ! Last index in z-direction

  REAL :: p  (nx,ny,nz)       ! Pressure (Pascal)
  REAL :: t  (nx,ny,nz)       ! Temperature (K)
  REAL :: qvsat(nx,ny,nz)     ! Saturation water vapor specific
                              ! humidity (kg/kg)
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$*  inline routine (f_qvsat)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=kbgn,kend
    DO j=jbgn,jend
      DO i=ibgn,iend
        qvsat(i,j,k) = f_qvsat( p(i,j,k), t(i,j,k) )
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE getqvs
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETMRS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getmrs( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend,             &,1
           p, t, mrsat )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the saturation mixing ratio using Teten's formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/09/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Array dimension in x-direction
!    ny       Array dimension in y-direction
!    nz       Array dimension in z-direction
!    ibgn     Starting index in x-direction
!    jbgn     Starting index in y-direction
!    kbgn     Starting index in z-direction
!    iend     Last index in x-direction
!    jend     Last index in y-direction
!    kend     Last index in z-direction
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!
!  OUTPUT:
!
!    mrsat    Saturation water vapor Mixing ratio (kg/kg).
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! Input array dimension

  INTEGER :: ibgn             ! Starting index in x-direction
  INTEGER :: jbgn             ! Starting index in y-direction
  INTEGER :: kbgn             ! Starting index in z-direction
  INTEGER :: iend             ! Last index in x-direction
  INTEGER :: jend             ! Last index in y-direction
  INTEGER :: kend             ! Last index in z-direction

  REAL :: p  (nx,ny,nz)       ! Pressure (Pascal)
  REAL :: t  (nx,ny,nz)       ! Temperature (K)
  REAL :: mrsat(nx,ny,nz)     ! Saturated water vapor mixing ratio (kg/kg)
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Function f_mrsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_mrsat

!fpp$ expand (f_mrsat)
!dir$ inline always f_mrsat
!*$*  inline routine (f_mrsat)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=kbgn,kend
    DO j=jbgn,jend
      DO i=ibgn,iend
        mrsat(i,j,k) = f_mrsat( p(i,j,k), t(i,j,k) )
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE getmrs
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETDEW                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getdew( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend,             & 2,1
           p, t, qv, tdew )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the dew point temperature (K) using enhanced Teten's
!  formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/09/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Array dimension in x-direction
!    ny       Array dimension in y-direction
!    nz       Array dimension in z-direction
!    ibgn     Starting index in x-direction
!    jbgn     Starting index in y-direction
!    kbgn     Starting index in z-direction
!    iend     Last index in x-direction
!    jend     Last index in y-direction
!    kend     Last index in z-direction
!
!    p        Pressure (Pascal)
!    t        Temperature (K)
!    qv       Specific humidity (kg/kg)
!
!  OUTPUT:
!
!    tdew     Dew point temperature (K)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! Input array dimension

  INTEGER :: ibgn             ! Starting index in x-direction
  INTEGER :: jbgn             ! Starting index in y-direction
  INTEGER :: kbgn             ! Starting index in z-direction
  INTEGER :: iend             ! Last index in x-direction
  INTEGER :: jend             ! Last index in y-direction
  INTEGER :: kend             ! Last index in z-direction

  REAL :: p  (nx,ny,nz)       ! Pressure (Pascal)
  REAL :: t  (nx,ny,nz)       ! Temperature (K)
  REAL :: qv (nx,ny,nz)       ! Specific humidity (kg/kg)
  REAL :: tdew(nx,ny,nz)      ! Dew point temperature (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Function f_mrsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_tdew

!fpp$ expand (f_tdew)
!dir$ inline always f_tdew
!*$*  inline routine (f_tdew)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=kbgn,kend
    DO j=jbgn,jend
      DO i=ibgn,iend
        tdew(i,j,k) = f_tdew( p(i,j,k), t(i,j,k), qv(i,j,k) )
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE getdew
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE PT2PTE                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE pt2pte( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend,             &,1
           p, pt, qv, pte )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the dew point temperature (K) using enhanced Teten's
!  formula.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/09/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Array dimension in x-direction
!    ny       Array dimension in y-direction
!    nz       Array dimension in z-direction
!    ibgn     Starting index in x-direction
!    jbgn     Starting index in y-direction
!    kbgn     Starting index in z-direction
!    iend     Last index in x-direction
!    jend     Last index in y-direction
!    kend     Last index in z-direction
!
!    p        Pressure (Pascal)
!    pt       Temperature (K)
!    qv       Specific humidity (kg/kg)
!
!  OUTPUT:
!
!    pte      Equivalent potential temperature (K)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! Input array dimension

  INTEGER :: ibgn             ! Starting index in x-direction
  INTEGER :: jbgn             ! Starting index in y-direction
  INTEGER :: kbgn             ! Starting index in z-direction
  INTEGER :: iend             ! Last index in x-direction
  INTEGER :: jend             ! Last index in y-direction
  INTEGER :: kend             ! Last index in z-direction

  REAL :: p  (nx,ny,nz)       ! Pressure (Pascal)
  REAL :: pt (nx,ny,nz)       ! Temperature (K)
  REAL :: qv (nx,ny,nz)       ! Specific humidity (kg/kg)
  REAL :: pte(nx,ny,nz)       ! Equivalent potential temperature (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Function f_pt2pte and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_pt2pte

!fpp$ expand (f_pt2pte)
!dir$ inline always f_pt2pte
!*$*  inline routine (f_pt2pte)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=kbgn,kend
    DO j=jbgn,jend
      DO i=ibgn,iend
        pte(i,j,k) = f_pt2pte( p(i,j,k), pt(i,j,k), qv(i,j,k) )
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE pt2pte
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_PCCL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_pccl(pm,p,t,td,mrbar,n),7
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  this function returns the pressure at the convective condensation
!  level given the appropriate sounding data.
!
!  the algorithm is decribed on p.17 of stipanuk, g.s.,1973:
!  "algorithms for generating a skew-t log p diagram and computing
!  selected meteorological quantities," atmospheric sciences labora-
!  tory, u.s. army electronics command, white sands missile range, new
!  mexico 88002.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    g.s. stipanuk     1973            original version.
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982
!
!
!  MODIFICATIONS:
!
!  01/30/1998 (Yuhe Liu)
!  Modified to use Kelvin and Pascal as units of temperature and
!  pressure
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    p       pressure (Pa). note that p(i).gt.p(i+1).
!    t       temperature (K)
!    td      dew point (K)
!    n       number of levels in the sounding and the dimension of
!            p, t and td
!    pm      pressure (Pa) at upper boundary of the layer for
!            computing the mean mixing ratio. p(1) is the lower
!            boundary.
!
!  OUTPUT:
!
!    mrbar   mean mixing ratio (kg/kg) in the layer bounded by
!            pressures p(1) at the bottom and pm at the top
!    f_pccl  pressure (Pa) at the convective condensation level
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: n

  REAL :: t(n),p(n),td(n)
  REAL :: pm

  REAL :: mrbar
  REAL :: f_pccl
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,l

  REAL :: pc,tq,del,x,a
  REAL :: mrsat1,mrsat2
!
!-----------------------------------------------------------------------
!
!  Function f_mrsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_mrsat
  REAL :: f_tmr

!fpp$ expand (f_mrsat)
!fpp$ expand (f_tmr)
!dir$ inline always f_mrsat, f_tmr
!*$*  inline routine (f_mrsat, f_tmr)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (pm /= p(1)) GO TO 5

  mrbar = f_mrsat( p(1),td(1) )
  pc= pm

  IF (ABS(  GO TO 25

  5   CONTINUE

  mrbar = 0.
  k = 0

  10   CONTINUE

  k = k+1
  IF (p(k) > pm) GO TO 10
  k = k-1
  j = k-1
  IF (j < 1) GO TO 20

!   compute the average mixing ratio....alog = natural log

  DO i= 1,j
    l = i+1
    mrsat1 = f_mrsat( p(i),td(i) )
    mrsat2 = f_mrsat( p(l),td(l) )
    mrbar = (mrsat1+mrsat2)*ALOG(  END DO

  20   CONTINUE
  l= k+1

!   estimate the dew point at pressure pm.

  mrsat1 = f_mrsat( p(k),td(k) )
  tq= td(k)+(td(l)-td(k))*(ALOG(pm/p(k)))/(ALOG(  mrsat2 = f_mrsat( pm,tq )
  mrbar=  mrbar+(mrsat1+mrsat2 )*ALOG(  mrbar=  mrbar/(2.*ALOG(

!   find level at which the mixing ratio line  mrbar crosses the
!   environmental temperature profile.

  25   CONTINUE

  DO j= 1,n
    i = n-j+1
    IF (p(i) < 300.) CYCLE

!   f_tmr = temperature (celsius) at pressure p (mb) along a mixing
!        ratio line given by  mrbar (g/kg)

    x = f_tmr(mrbar,p(i))-t(i)
    IF (x <= 0.) GO TO 35
  END DO

  f_pccl = 0.0

  RETURN

!  set up bisection routine

  35   l = i
  i = i+1
  del = p(l)-p(i)
  pc = p(i)+.5*del
  a = (t(i)-t(l))/ALOG(  DO j = 1,10
    del = del/2.
    x = f_tmr(mrbar,pc)-t(l)-a*(ALOG(
!   the sign function replaces the sign of the first argument
!   with that of the second.

    pc = pc+SIGN(del,x)
  END DO

  45   f_pccl = pc

  RETURN
  END FUNCTION f_pccl
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_CT                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_ct(mrbar,pc,ps),1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This function returns the convective temperature ct (K)
!  given the mean mixing ratio mrbar (kg/kg) in the surface layer,
!  the pressure pc (Pa) at the convective condensation level (ccl)
!  and the surface pressure ps (Pa).
!
!-----------------------------------------------------------------------
!
!  AUTHOR: G.S. Stipanuk     1973            original version.
!
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982
!
!  MODIFICATION HISTORY:
!
!  01/30/1998 (Yuhe Liu)
!  Modified for ARPS to use the units of Kelvin, Pascal and kg/kg for
!  temperature, pressure and mixing ratio respectively
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    mrbar    Mean mixing ratio (kg/kg)
!    pc       Pressure at CCL (Pa)
!    ps       Pressure at surface (Pa)
!
!  OUTPUT:
!
!    f_ct     Convective temperature (K)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: mrbar
  REAL :: pc
  REAL :: ps

  REAL :: f_ct

  REAL :: tc
!
!-----------------------------------------------------------------------
!
!  Include file
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_tmr and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_tmr

!fpp$ expand (f_tmr)
!dir$ inline always f_tmr
!*$*  inline routine (f_tmr)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  compute the temperature (K) at the ccl.

  tc= f_tmr(mrbar,pc)
!
!-----------------------------------------------------------------------
!
!  compute the potential temperature (K) at ccl, i.e., the dry
!  adiabat through the ccl.
!
!  ptc = tc*((100000./pc)**rddcp)
!
!  and then compute the surface temperature on the same dry adiabat.
!
!  f_ct = ptc*((ps/100000.)**rddcp)
!
!-----------------------------------------------------------------------
!
  f_ct = tc*((ps/pc)**rddcp)

  RETURN
  END FUNCTION f_ct
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION F_TMR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION f_tmr(p,mr) 3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This function returns the temperature (Kelvin) on a mixing
!  ratio line mr (kg/kg) at pressure p (Pa). The formula is given in
!  table 1 on page 7 of Stipanuk (1973).
!
!-----------------------------------------------------------------------
!
!  AUTHOR: G.S. Stipanuk, 1973            original version.
!
!    reference stipanuk paper entitled:
!         "algorithms for generating a skew-t, log p
!         diagram and computing selected meteorological
!         quantities."
!         atmospheric sciences laboratory
!         u.s. army electronics command
!         white sands missile range, new mexico 88002
!         33 pages
!    baker, schlatter  17-may-1982
!
!  MODIFICATIONS:
!
!  01/30/1998 (Yuhe Liu)
!  Modified to use the units of Kelvin, Pascal and kg/kg for
!  temperature, pressure and mixing ratio respectively.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    p        Pressure (Pa)
!    mr       Mixing ratio (kg/kg)
!
!  OUTPUT:
!
!    f_tmr    Temperature (K)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: p              ! Pressure (Pa)
  REAL :: mr             ! Mixing ratio (kg/kg)

  REAL :: f_tmr          ! Temperature (K)
!
!-----------------------------------------------------------------------
!
!  Local variables
!
!-----------------------------------------------------------------------
!
  REAL :: x, pmb, mrg

  REAL :: c1,c2,c3,c4,c5,c6

  DATA c1/.0498646455/
  DATA c2/2.4082965/
  DATA c3/7.07475/
  DATA c4/38.9114/
  DATA c5/.0915/
  DATA c6/1.2035/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  pmb = p/100.
  mrg = mr*1000.

  x= ALOG10(mrg*pmb/(622.+mrg))
  f_tmr= 10.**(c1*x+c2)-c3+c4*((10.**(c5*x)-c6)**2.)

  RETURN
  END FUNCTION f_tmr