real function Cefp( qv, dl )
! 水蒸気混合比から有効定圧比熱を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル

  Cefp=(Cpd+qv*Cpv)/(1.0+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'Cefp', Cefp, 'J K-1 kg-1' )
  end if

  return
end function

real function Cefv( qv, dl )
! 水蒸気混合比から有効定積比熱を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル

  Cefv=(Cvd+qv*Cvv)/(1.0+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'Cefv', Cefv, 'J K-1 kg-1' )
  end if

  return
end function

real function Cl( T, dl )
! 液水の比熱を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: c1=4190.0, c2=4770.0
  real :: val_rat, tref

  val_rat=(c1-c2)/40.0
  tref=273.15-40.0

  if(T>=273.15)then
     Cl=c1
  else
     Cl=c2+val_rat*(T-tref)
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'Cl', Cl, 'J K-1 kg-1' )
  end if

  return
end function

real function DSE_Emanuel( T, qv, z, dl )
! Emanuel (1994) の (4.5.24) で定義される乾燥静的エネルギーを計算する.
  use Thermo_Const
  use Phys_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in) :: z  ! 地表面高度 [m]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n

  DSE_Emanuel=(Cpd+qv*Cpv)*T+(1.0+qv)*g*z

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'DSE_Emanuel', DSE_Emanuel, 'J kg-1' )
  end if

  return
end function

real function LH(T, dl)  ! 温度 T における潜熱の計算
  ! 本関数は, 潜熱を一定とした場合の結果とそれほど大きな差は生まれない.
  ! 本計算式は, 潜熱の温度変化に関する方程式「キルヒホッフの式」を用いた.
  ! また, その際に必要な液相の水の定圧比熱および水蒸気圧の定圧比熱は
  ! 温度依存性がないものと仮定し, それぞれ
  ! $C_l=4190,\; C_{pv}=1870$という値を用いて導出した関係式である.
  ! よって, 用いる比熱の値やその温度依存性を考慮すると係数が少し変化する可能性.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル

  LH=LH0-2.32e3*(T-t0)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'LH', LH, 'J kg-1' )
  end if

  return
end function

real function MSE_Emanuel( T, qv, z, qo, dl )
! Emanuel (1994) の (4.5.23) で定義される湿潤静的エネルギーを計算する.
  use Thermo_Const
  use Phys_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in) :: z  ! 地表面高度 [m]
  real, intent(in), optional :: qo(:)  ! 水蒸気以外の水物質混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq

  if(present(qo))then
     n=size(qo)
     tmpq=qv
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=qv
  end if

  MSE_Emanuel=(Cpd+tmpq*Cl( T ))*T+LH( T )*qv+(1.0+tmpq)*g*z

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'MSE_Emanuel', MSE_Emanuel, 'J kg-1' )
  end if

  return
end function

real function RHTP_2_qv(RH,T,P, dl)  ! 相対湿度と温度から混合比を計算する
  ! RHT_2_e から水蒸気圧を計算し, eP_2_qv から混合比を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: RH  ! 相対湿度 [%]
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: e

  e=RHT_2_e(RH,T)
  RHTP_2_qv=eP_2_qv(e,P)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'RHTP_2_qv', RHTP_2_qv, 'kg kg-1' )
  end if

  return
end function

real function RHT_2_e(RH,T, dl)  ! 相対湿度と温度から水蒸気圧を計算する
  ! $RH=(e/es)\times 100$ という定義から計算.
  use Thermo_Const
  implicit none
  real, intent(in) :: RH  ! 相対湿度 [%]
  real, intent(in) :: T  ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: es

  es=es_Bolton(T)
  RHT_2_e=RH*es*1.0e-2

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'RHT_2_e', RHT_2_e, 'Pa' )
  end if

  return
end function

real function Reff( qv, dl )
! 水蒸気混合比から有効気体定数を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル

  Reff=(Rd+qv*Rv)/(1.0+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'Reff', Reff, 'J K-1 kg-1' )
  end if

  return
end function

real function TP_2_qvs(T,P, dl)  ! 温度と全圧から飽和混合比を計算する
  ! ここでは, es_Bolton を用いて飽和水蒸気圧を計算した後,
  ! eP_2_qv を用いて混合比に変換することで飽和混合比を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 大気の全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps
  real :: es

  eps=Rd/Rv
  es=es_Bolton(T)
  TP_2_qvs=eps*es/(P-es)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TP_2_qvs', TP_2_qvs, 'kg kg-1' )
  end if

  return
end function

real function TP_2_rho( T, P, dl )
! 乾燥大気の状態方程式から, 温度と気圧を与えて密度を得る.
  use Thermo_Const
  implicit none
  real, intent(in) :: T    ! 大気の温度 [K]
  real, intent(in) :: P    ! 大気の圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル

  TP_2_rho=p/(Rd*T)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TP_2_rho', TP_2_rho, 'kg m-3' )
  end if

  return
end function

real function TTd_2_RH_Bolton( T, Td, dl )
  ! 温度と露点温度から相対湿度を計算する.
  ! Bolton の式を用いて水蒸気圧・飽和水蒸気圧を計算する. 
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: Td ! 露点温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: e,es

  e=es_Bolton(Td)
  es=es_Bolton(T)

  TTd_2_RH_Bolton=(e/es)*100.0

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TTd_2_RH_Bolton', TTd_2_RH_Bolton, '%' )
  end if
  
  return
end function

real function TTd_2_RH_tetens( T, Td, dl )
  ! 温度と露点温度から相対湿度を計算する.
  ! tetens の式を用いて水蒸気圧・飽和水蒸気圧を計算する. 
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: Td ! 露点温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: e,es

  e=tetens(Td)
  es=tetens(T)

  TTd_2_RH_tetens=(e/es)*100.0

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TTd_2_RH_tetens', TTd_2_RH_tetens, '%' )
  end if
  
  return
end function

real function Tq_2_Trho( T, qv, qo, dl )
! 温度と水物質から密度温度を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in), optional :: qo(:)  ! 水蒸気以外の水凝結物質混合比 [kg/kg]
                   ! 凝結物質のカテゴリはいくらでも構わない.
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq

  if(present(qo))then
     n=size(qo)
     tmpq=qv
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=qv
  end if

  Tq_2_Trho=T*(1.0+qv*Rv/Rd)/(1.0+tmpq)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'Tq_2_Trho', Tq_2_Trho, 'K' )
  end if

  return
end function

real function TqvP_2_TLCL(T,qv,P, dl)  !! 温度と混合比と全圧から T_LCL を計算する
  ! 混合比から水蒸気圧を求め, そこから T_LCL を計算する
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 混合比 [kg / kg]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: coe=2840.0, a=3.5, b=4.805, c=55.0
  real :: e

  e=qvP_2_e(qv,P)
  e=e*1.0e-2
  if(e>0.0)then
     TqvP_2_TLCL=coe/(a*log(T)-log(e)-b)+55.0
  else
     TqvP_2_TLCL=coe/(a*log(T)-b)+55.0  ! true ??
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TqvP_2_TLCL', TqvP_2_TLCL, 'K' )
  end if

  return
end function

real function TqvP_2_thetae(T,qv,P, dl)  ! 温度, 混合比, 全圧から相当温位を計算する.
  ! T_LCL を用いるので, そのための関数を使用する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 混合比 [kg / kg]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: T_LCL, kappa, theta_d, e, qvs
  real, parameter :: a=0.2854, b=0.28, c=3376.0, d=0.81

  kappa=Rd/Cpd
  e=qvP_2_e(qv,P)
  T_LCL=TqvP_2_TLCL(T,qv,P)
  theta_d=T*(p0/(P-e))**kappa
  qvs=TP_2_qvs(T_LCL,P)
  TqvP_2_thetae=theta_d*exp(LH(T_LCL)*qvs/(Cpd*T_LCL))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TqvP_2_thetae', TqvP_2_thetae, 'K' )
  end if

  return
end function

real function TqvP_2_thetaes(T,P, dl)  ! 温度, 混合比, 全圧から飽和相当温位を計算する.
  ! T_LCL を用いるので, そのための関数を使用する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
!  real, intent(in) :: qv  ! 混合比 [kg / kg]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa, theta_d, qvs !,e
  real, parameter :: a=0.2854, b=0.28, c=3376.0, d=0.81

  kappa=Rd/Cpd
!  e=qvP_2_e(qv,P)
!  theta_d=T*(p0/(P-e))**kappa
  theta_d=T*(p0/P)**kappa
  qvs=TP_2_qvs(T,P)
  TqvP_2_thetaes=theta_d*exp(LH(T)*qvs/(Cpd*T))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TqvP_2_thetaes', TqvP_2_thetaes, 'K' )
  end if

  return
end function

real function TqvP_2_thetav( T, qv, P, dl )  ! 温度, 水蒸気混合比, 圧力から仮温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg / kg]
  real, intent(in) :: T   ! 温度 [K]
  real, intent(in) :: P   ! 圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa, Tv

  kappa=Rd/cpd
  Tv=qvT_2_Tv(qv,T)
  TqvP_2_thetav=theta_dry(Tv,P)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TqvP_2_thetav', TqvP_2_thetav, 'K' )
  end if

  return
end function

real function TthetavP_2_qv( T, thetav, P, dl )  ! 温度, 仮温位, 圧力から水蒸気混合比を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T       ! 温度 [K]
  real, intent(in) :: thetav  ! 仮温位 [K]
  real, intent(in) :: P       ! 圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: Tv, eps, exn

  eps=Rd/Rv

  exn=exner_func_dry( P )

  if(T/=exn*thetav)then
     TthetavP_2_qv=eps*((T-exn*thetav)/(eps*exn*thetav-T))
  else
     TthetavP_2_qv=0.0
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TthetavP_2_qv', TthetavP_2_qv, 'kg kg-1' )
  end if

  return
end function

real function TvT_2_qv( Tv, T, dl)  ! 仮温度と温度から水蒸気混合比を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: Tv  ! 仮温度 [K]
  real, intent(in) :: T   ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps

  eps=Rd/Rv

  if(Tv/=T)then
     TvT_2_qv=eps*((T-Tv)/(eps*Tv-T))
  else
     TvT_2_qv=0.0
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'TvT_2_qv', TvT_2_qv, 'kg kg-1' )
  end if

  return
end function

real function WSE_Emanuel( T, ql, z, qo, dl )
! Emanuel (1994) の (4.5.25) で定義される液水静的エネルギーを計算する.
  use Thermo_Const
  use Phys_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: ql  ! 液水混合比 [kg/kg]
  real, intent(in) :: z  ! 地表面高度 [m]
  real, intent(in), optional :: qo(:)  ! 液水以外の水物質混合比（水蒸気含む） [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq

  if(present(qo))then
     n=size(qo)
     tmpq=ql
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=ql
  end if

  WSE_Emanuel=(Cpd+tmpq*Cpv)*T-LH( T )*ql+(1.0+tmpq)*g*z

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'WSE_Emanuel', WSE_Emanuel, 'J kg-1' )
  end if

  return
end function

real function eP_2_qv(e,P, dl)  ! 水蒸気圧と全圧から混合比を計算する
  use Thermo_Const
  implicit none
  real, intent(in) :: e  ! 水蒸気圧 [Pa]
  real, intent(in) :: P  ! 大気の全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps

  eps=Rd/Rv
  eP_2_qv=eps*e/(P-e)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'eP_2_qv', eP_2_qv, 'kg kg-1' )
  end if

  return
end function

real function eT_2_RH(e,T, dl)  ! 水蒸気圧と温度から相対湿度を計算する
  ! $RH=(e/es)\times 100$ という定義から計算.
  use Thermo_Const
  implicit none
  real, intent(in) :: e  ! 水蒸気圧 [Pa]
  real, intent(in) :: T  ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: es

  es=es_Bolton(T)
  eT_2_RH=100.0*e/es

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'eT_2_RH', eT_2_RH, '%' )
  end if

  return
end function

real function es_Bolton(T, dl)  ! Bolton(1980) の手法を用いて飽和水蒸気圧を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=17.67, c=29.65

  es_Bolton=e0*exp(a*((T-t0)/(T-c)))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'es_Bolton', es_Bolton, 'Pa' )
  end if

  return
end function

real function es_TD(e, dl)  ! es_Bolton を用いて飽和水蒸気圧の計算式の逆算から
! 露点温度を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: e  ! 大気の水蒸気圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=17.67, c=29.65

  if(e>0.0)then
     es_TD=c+(a*(t0-c))/(a-log(e/e0))
  else
     es_TD=t0  ! true ??
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'es_TD', es_TD, 'Pa' )
  end if

  return
end function

real function esi_Emanuel( T, dl )
! Emanuel (1994) の (4.4.15) から氷飽和蒸気圧を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=23.33086, b=6111.72784, c=0.15215

  esi_Emanuel=a-b/T+c*log(T)
  esi_Emanuel=exp(esi_Emanuel)*100.0

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'esi_Emanuel', esi_Emanuel, 'Pa' )
  end if

  return
end function

real function exner_func_dry(P, dl)  ! 乾燥大気についてのエクスナー関数を計算する
  use Thermo_Const
  implicit none
  real, intent(in) :: P  ! 圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa

  kappa=Rd/Cpd
  exner_func_dry=(P/p0)**kappa

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'exner_func_dry', exner_func_dry, '1' )
  end if

  return
end function

real function get_gamma_d( dl )  ! 乾燥断熱減率を呼ぶ関数(この関数必要か？)
  use Thermo_Const
  use Phys_Const
  implicit none
  integer, intent(in), optional :: dl  ! デバッグレベル

  get_gamma_d=-g/Cpd

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'get_gamma_d', get_gamma_d, 'K m-1' )
  end if

  return
end function

real function goff_gratch(T, dl)  ! goff-gratch の式を用いて飽和水蒸気圧を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=-7.90298, b=5.02808, c=-1.3816e-7, d=8.1328e-3
  real, parameter :: pa=11.344, pb=-3.49149
  real, parameter :: tst=373.15
  real, parameter :: est=1.01325e5
  real :: term

  term=a*(tst/T-1.0)+b*log10(tst/T)+c*(10.0**(pa*(1.0-T/tst))-1.0)+d*(10.0**(pb*(tst/T-1.0))-1.0)
  goff_gratch=est*10.0**(term)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'goff_gratch', goff_gratch, 'Pa' )
  end if

  return
end function

real function goff_gratch_i(T, dl)  ! goff-gratch の式を用いて飽和水蒸気圧 (氷飽和) を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=-9.09718, b=-3.56654, c=0.876793
  real, parameter :: est=6.1173e2
  real :: term

  term=a*(ti0/T-1.0)+b*log10(ti0/T)+c*(1.0-t/ti0)
  goff_gratch_i=ei0*10.0**(term)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'goff_gratch_i', goff_gratch_i, 'Pa' )
  end if

  return
end function

real function hypsometric_form(p,z,T,z_t, dl)  ! 高度と気圧を与えてある高度の気圧を求める.
!  気圧の更正には対流圏における標準大気の温度減率である 6.5 [K/km] で計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: p  ! 基準となる高度における圧力 [Pa]
  real, intent(in) :: z  ! 基準となる高度 [m]
  real, intent(in) :: T  ! 基準点での温度 [K]
  real, intent(in), optional :: z_t  ! 求める高度 [m] : デフォルトでは 0 m.
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: gam = 6.5e-3  ! 標準大気の温度減率 [K/m]
  real :: z_calc, p_tmp

!write(*,*) "hypsometric, g=", g

  if(present(z_t))then
     z_calc=z_t
  else
     z_calc=0.0
  end if

  p_tmp=p*((T+gam*z)/(T+gam*z_calc))**(g/(gam*Rd))
  hypsometric_form=p_tmp

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'hypsometric_form', hypsometric_form, 'Pa' )
  end if

  return
end function

real function liquid_enthal( T, ql, qo, dl )
! 温度と水物質から液水エンタルピーを計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: ql  ! 液体水混合比 [kg/kg]
  real, intent(in), optional :: qo(:)  ! 液体水以外の水物質混合比 [kg/kg]
                   ! 水物質のカテゴリはいくらでも構わない（水蒸気も含む）.
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq

  if(present(qo))then
     n=size(qo)
     tmpq=ql
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=ql
  end if

  liquid_enthal=(Cpd+tmpq*Cpv)*T-LH( T )*ql

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'liquid_enthal', liquid_enthal, 'J kg-1' )
  end if

  return
end function

real function moist_enthal( T, qv, qo, dl )
! 温度と水物質から湿潤エンタルピーを計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in), optional :: qo(:)  ! 水蒸気以外の水凝結物質混合比 [kg/kg]
                   ! 凝結物質のカテゴリはいくらでも構わない.
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq

  if(present(qo))then
     n=size(qo)
     tmpq=qv
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=qv
  end if

  moist_enthal=(Cpd+tmpq*Cl( T ))*T+LH( T )*qv

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'moist_enthal', moist_enthal, 'J kg-1' )
  end if

  return
end function

real function moist_laps_temp( p_stan, T_stan, p, dl )
! 偽断熱過程における湿潤断熱減率に従う, 高度 p における気温を求める関数.
! 気温はパーセルの気温.
  use Thermo_Const
  implicit none
  real, intent(in) :: p_stan  ! 基準気圧座標 [Pa]
  real, intent(in) :: T_stan  ! 基準温度 [K]
  real, intent(in) :: p  ! 求める高度における気圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: seqpt, t_temp, seqpt_temp, dt, seqpt_temp2
  integer, parameter :: ilim=20
  real, parameter :: limit=2.0e-4

  seqpt=thetaes_Bolton( T_stan, p_stan )
  seqpt_temp=thetaes_Bolton( t_stan, p )
  if(p>p_stan)then
     dt=1.0
  else
     dt=-1.0
  end if
  t_temp=t_stan+dt
  seqpt_temp2=thetaes_Bolton( t_temp, p )

  do while ( (seqpt_temp2-seqpt)*(seqpt_temp-seqpt) >0.0 )
     t_temp=t_temp+dt
     seqpt_temp=seqpt_temp2
     seqpt_temp2=thetaes_Bolton( t_temp, p )
  end do

  seqpt_temp=seqpt_temp2
  dt=1.0

  do while ( abs(seqpt_temp-seqpt) > limit )

     dt=0.5*dt

     if( (seqpt_temp-seqpt)>0.0 )then
        t_temp=t_temp-dt
     else
        t_temp=t_temp+dt
     end if

     seqpt_temp=thetaes_Bolton( t_temp, p )

  end do

  moist_laps_temp=t_temp

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'moist_laps_temp', moist_laps_temp, 'K' )
  end if

  return
end function

real function qvP_2_e(qv,P, dl)  ! 混合比と全圧から水蒸気圧を計算する
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 混合比 [kg / kg]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps

  eps=Rd/Rv
  qvP_2_e=P*qv/(eps+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'qvP_2_e', qvP_2_e, 'Pa' )
  end if

  return
end function

real function qvTP_2_RH(qv,T,P, dl)  ! 混合比と温度から相対湿度を計算する.
  ! qvP_2_e から水蒸気圧を計算し, 相対湿度の定義を用いる.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 相対湿度 [kg / kg]
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: e

  e=qvP_2_e(qv,P)
  qvTP_2_RH=eT_2_RH(e,T)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'qvTP_2_RH', qvTP_2_RH, '%' )
  end if

  return
end function

real function qvT_2_Tv( qv, T, dl)  ! 温度と水蒸気混合比から仮温度を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg / kg]
  real, intent(in) :: T   ! 温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps

  eps=Rd/Rv
  qvT_2_Tv=T*(1.0+qv/eps)/(1.0+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'qvT_2_Tv', qvT_2_Tv, 'K' )
  end if

  return
end function

real function qvTv_2_T( qv, Tv, dl )  ! 水蒸気混合比と仮温度から温度を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg / kg]
  real, intent(in) :: Tv  ! 仮温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps

  eps=Rd/Rv
  qvTv_2_T=Tv*(1.0+qv)/(1.0+qv/eps)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'qvTv_2_T', qvTv_2_T, 'K' )
  end if

  return
end function

real function qv_2_sh( qv, dl )
! 混合比から比湿を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル

  qv_2_sh=qv/(1.0+qv)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'qv_2_sh', qv_2_sh, 'kg kg-1' )
  end if

  return
end function

real function rhoP_2_T( rho, P, dl )
! 乾燥大気の状態方程式から, 密度と気圧を与えて温度を得る.
  use Thermo_Const
  implicit none
  real, intent(in) :: rho  ! 大気の密度 [kg/m^3]
  real, intent(in) :: P    ! 大気の圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル

  rhoP_2_T=P/(Rd*rho)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'rhoP_2_T', rhoP_2_T, 'K' )
  end if

  return
end function

real function rhoT_2_P( rho, T, dl )
! 乾燥大気の状態方程式から, 密度と温度を与えて気圧を得る.
  use Thermo_Const
  implicit none
  real, intent(in) :: rho  ! 大気の密度 [kg/m^3]
  real, intent(in) :: T    ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル

  rhoT_2_P=rho*Rd*T

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'rhoT_2_P', rhoT_2_P, 'Pa' )
  end if

  return
end function

real function rho_ocean( pres, temp, sal, dl )
! UNESCO (1981) の状態方程式を用いて, 絶対密度を計算する.
! ここでの計算式には以下のようなコード作成上の処理が含まれていることに注意する.
!    (a**x)*(b**y)=exp(x*ln(a)+y*ln(b))
! なぜなら, 
!    (a**x)*(b**y)=z
! => ln{(a**x)*(b**y)}=ln(z)
! => x*ln(a)+y*ln(b)=ln(z)
! => exp(x*ln(a)+y*ln(b))=z.
! とくに, c*(b**y) という形の場合 (係数 x 変数べき乗) :
! c*(b**y)=exp(ln(c)+y*ln(b)) で表される.
  implicit none
  real, intent(in) :: pres    ! 圧力 [Bar]
  real, intent(in) :: temp    ! 温度 [degC]
  real, intent(in) :: sal     ! 塩分濃度 [PSU]
  integer, intent(in), optional :: dl  ! デバッグレベル

  real :: rhoz0, rhowc, Kw, Kz0, Kz, tmptemp

  if(temp>0.0)then
     tmptemp=temp
     Kw=19652.21+148.4206*tmptemp-2.327105*tmptemp*tmptemp +exp(log(1.360477e-2)+3.0*log(tmptemp)) -exp(log(5.155288e-5)+4.0*log(tmptemp))
     Kz0=Kw+sal*(54.6746-0.603459*tmptemp +exp(log(1.09987e-2)+2.0*log(tmptemp)) -exp(log(6.167e-5)+3.0*log(tmptemp))) +(sal**1.5)*(7.944e-2+1.6483e-2*tmptemp -exp(log(5.3009e-4)+2.0*log(tmptemp)))

     Kz=Kz0+pres*(3.239908+1.43713e-3*tmptemp +exp(log(1.16092e-4)+2.0*log(tmptemp)) -exp(log(5.77905e-7)+3.0*log(tmptemp))) +pres*sal*(2.2838e-3-1.0981e-5*tmptemp -exp(log(1.6078e-6)+2.0*log(tmptemp))) +pres*(sal**1.5)*1.91075e-4 +(pres**2)*(8.50935e-5-6.12293e-6*tmptemp +exp(log(5.2787e-8)+2.0*log(tmptemp))) +(pres**2)*sal*(-9.9348e-7+2.0816e-8*tmptemp +exp(log(9.1697e-10)+2.0*log(tmptemp)))

     rhowc=999.842594+6.793952e-2*tmptemp -exp(log(9.09529e-3)+2.0*log(tmptemp)) +exp(log(1.001685e-4)+3.0*log(tmptemp)) -exp(log(1.120083e-6)+4.0*log(tmptemp)) +exp(log(6.536332e-9)+5.0*log(tmptemp))

     rhoz0=rhowc +sal*(0.824493-4.0899e-3*tmptemp +exp(log(7.6438e-5)+2.0*log(tmptemp)) -exp(log(8.2467e-7)+3.0*log(tmptemp)) +exp(log(5.3875e-9)+4.0*log(tmptemp))) +(sal**1.5)*(-5.72466e-3+1.0227e-4*tmptemp -exp(log(1.6546e-6)+2.0*log(tmptemp))) +(sal**2)*4.8314e-4

  else if(temp<0.0)then

     tmptemp=-temp
     Kw=19652.21-148.4206*tmptemp-2.327105*tmptemp*tmptemp -exp(log(1.360477e-2)+3.0*log(tmptemp)) -exp(log(5.155288e-5)+4.0*log(tmptemp))

     Kz0=Kw+sal*(54.6746+0.603459*tmptemp +exp(log(1.09987e-2)+2.0*log(tmptemp)) +exp(log(6.167e-5)+3.0*log(tmptemp))) +(sal**1.5)*(7.944e-2-1.6483e-2*tmptemp -exp(log(5.3009e-4)+2.0*log(tmptemp)))

     Kz=Kz0+pres*(3.239908-1.43713e-3*tmptemp +exp(log(1.16092e-4)+2.0*log(tmptemp)) +exp(log(5.77905e-7)+3.0*log(tmptemp))) +pres*sal*(2.2838e-3+1.0981e-5*tmptemp -exp(log(1.6078e-6)+2.0*log(tmptemp))) +pres*(sal**1.5)*1.91075e-4 +(pres**2)*(8.50935e-5+6.12293e-6*tmptemp +exp(log(5.2787e-8)+2.0*log(tmptemp))) +(pres**2)*sal*(-9.9348e-7-2.0816e-8*tmptemp +exp(log(9.1697e-10)+2.0*log(tmptemp)))

     rhowc=999.842594-6.793952e-2*tmptemp -exp(log(9.09529e-3)+2.0*log(tmptemp)) -exp(log(1.001685e-4)+3.0*log(tmptemp)) -exp(log(1.120083e-6)+4.0*log(tmptemp)) -exp(log(6.536332e-9)+5.0*log(tmptemp))

     rhoz0=rhowc +sal*(0.824493+4.0899e-3*tmptemp +exp(log(7.6438e-5)+2.0*log(tmptemp)) +exp(log(8.2467e-7)+3.0*log(tmptemp)) +exp(log(5.3875e-9)+4.0*log(tmptemp))) +(sal**1.5)*(-5.72466e-3-1.0227e-4*tmptemp -exp(log(1.6546e-6)+2.0*log(tmptemp))) +(sal**2)*4.8314e-4

  else

     Kw=19652.21
     Kz0=Kw+sal*(54.6746 +(sal**1.5)*7.944e-2)
     Kz=Kz0+pres*3.239908 +pres*sal*2.2838e-3 +pres*(sal**1.5)*1.91075e-4 +(pres**2)*8.50935e-5 +(pres**2)*sal*(-9.9348e-7)

     rhowc=999.842594

     rhoz0=rhowc +sal*0.824493 -5.72466e-3*(sal**1.5) +4.8314e-4*(sal**2)

  end if

  rho_ocean=rhoz0/(1.0-pres/Kz)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'rho_ocean', rho_ocean, 'kg m-3' )
  end if

  return
end function rho_ocean

real function sh_2_qv( sh, dl )
! 比湿から水蒸気混合比を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: sh  ! 比湿 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル

  sh_2_qv=sh/(1.0-sh)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'sh_2_qv', sh_2_qv, 'kg kg-1' )
  end if

  return
end function

real function tetens( T, dl )  ! テテンの実験式を用いて飽和水蒸気圧を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 大気の温度 [K]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: a=7.5, b=237.7, c=9.5, d=265.5

  if(t<=t0)then
     tetens=e0*10.0**(c*(t-t0)/(t-t0+d))
  else
     tetens=e0*10.0**(a*(t-t0)/(t-t0+b))
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'tetens', tetens, 'Pa' )
  end if

  return
end function

real function thetaP_2_T(theta,P, dl)  ! 温位, 圧力から温度を計算する(乾燥大気として計算)
  use Thermo_Const
  implicit none
  real, intent(in) :: theta  ! 温位 [K]
  real, intent(in) :: P  ! 湿潤大気の全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa

  kappa=Rd/Cpd

  thetaP_2_T=theta*(P/p0)**kappa

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetaP_2_T', thetaP_2_T, 'K' )
  end if

  return
end function

real function thetaT_2_P(theta,T, dl)  ! 温位, 温度から圧力を計算する(乾燥大気として計算)
  use Thermo_Const
  implicit none
  real, intent(in) :: theta  ! 温位 [K]
  real, intent(in) :: T  ! 温度 [T]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa

  kappa=Cpd/Rd

  thetaT_2_P=p0*(T/theta)**kappa

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetaT_2_P', thetaT_2_P, 'Pa' )
  end if

  return
end function

real function theta_dry(T,P, dl)  ! 乾燥大気における温位を計算する
  ! ただし, 湿潤大気においても, 観測される全圧を P として計算することができ
  ! その結果は別関数 theta_moist の結果とそれほど変わらない.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 乾燥大気の気圧(もしくは, 湿潤大気の全圧) [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa

  kappa=Rd/Cpd
  theta_dry=T*(p0/P)**kappa

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'theta_dry', theta_dry, 'K' )
  end if

  return
end function

real function theta_moist(T,P,qv, dl)  ! 湿潤大気における温位を計算する
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 湿潤大気の全圧 [Pa]
  real, intent(in) :: qv  ! 混合比 [kg / kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: eps, kappa, CR

  eps=Rd/Rv
  kappa=Rd/Cpd
  CR=Cpv/Cpd

  kappa=kappa*((1.0+qv*CR)/(1.0+qv/eps))  ! kappa の値が上から上書きされていることに注意
  theta_moist=T*(p0/P)**kappa

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'theta_moist', theta_moist, 'K' )
  end if

  return
end function

real function thetae_Bolton(T,qv,P, dl)
  ! Bolton(1980) による手法を用いて相当温位を計算する.
  ! この相当温位は偽断熱過程での相当温位である.
  ! T_LCL を用いるので, そのための関数を使用する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 混合比 [kg / kg]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: T_LCL, qvs
  real, parameter :: a=0.2854, b=0.28, c=3376.0, d=0.81

  T_LCL=TqvP_2_TLCL(T,qv,P)
  qvs=TP_2_qvs(T_LCL,P)
  thetae_Bolton=T*((p0/P)**(a*(1.0-b*qvs))) *exp((c/T_LCL-2.54)*qvs*(1.0+d*qvs))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetae_Bolton', thetae_Bolton, 'K' )
  end if

  return
end function

real function thetae_Emanuel( T, qv, pres, qo, dl )
! Emanuel (1994) の (4.5.11) で定義される相当温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in) :: pres  ! 圧力 [Pa]
  real, intent(in), optional :: qo(:)  ! 水蒸気以外の水物質混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq, calH, pdry, pow1, pow2, pow3, coe1, coe2, coe3, Cltmp, LHtmp

  if(present(qo))then
     n=size(qo)
     tmpq=qv
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=qv
  end if

  pdry=pres-qvP_2_e( qv, pres )
  calH=qvTP_2_RH( qv, T, pres )*0.01
  Cltmp=Cl( T )
  LHtmp=LH( T )

  pow1=Rd/(Cpd+Cltmp*tmpq)
  pow2=qv*Rv/(Cpd+Cltmp*tmpq)
  pow3=LHtmp*qv/((Cpd+Cltmp*tmpq)*T)
  coe1=(p0/pdry)**pow1
  coe2=calH**pow2
  coe3=exp(pow3)

  thetae_Emanuel=T*coe1*coe2*coe3

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetae_Emanuel', thetae_Emanuel, 'K' )
  end if

  return
end function

real function thetaes_Bolton(T,P, dl)  ! Bolton(1980) による手法を用いて飽和相当温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: P  ! 全圧 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: qvs
  real, parameter :: a=0.2854, b=0.28, c=3376.0, d=0.81

  qvs=TP_2_qvs(T,P)
  thetaes_Bolton=T*((p0/P)**(a*(1.0-b*qvs))) *exp((c/T-2.54)*qvs*(1.0+d*qvs))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetaes_Bolton', thetaes_Bolton, 'K' )
  end if

  return
end function

real function thetal_Emanuel( T, ql, pres, qo, dl )
! Emanuel (1994) の (4.5.15) で定義される液水温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: ql  ! 液水混合比 [kg/kg]
  real, intent(in) :: pres  ! 圧力 [Pa]
  real, intent(in), optional :: qo(:)  ! 液水以外の水物質混合比（水蒸気含む） [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq, pow1, pow2, pow3, coe1, coe2, coe3, epsi

  epsi=Rd/Rv

  if(present(qo))then
     n=size(qo)
     tmpq=ql
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=ql
  end if

  pow1=(Rd+tmpq*Rv)/(Cpd+tmpq*Cpv)
  pow2=(tmpq*Rv)/(Cpd+tmpq*Cpv)
  pow3=-LH( T )*ql/((Cpd+tmpq*Cpv)*T)
  coe1=((p0/pres)*(1.0-ql/(epsi+tmpq)))**pow1
  coe2=(tmpq/(tmpq-ql))**pow2
  coe3=exp(pow3)

  thetal_Emanuel=T*coe1*coe2*coe3

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetal_Emanuel', thetal_Emanuel, 'K' )
  end if

  return
end function

real function thetalv_Emanuel( T, ql, qv, pres, qo, dl )
! Emanuel (1994) の (4.5.18) で定義される液水仮温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: ql  ! 液水混合比 [kg/kg]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in) :: pres  ! 圧力 [Pa]
  real, intent(in), optional :: qo(:)  ! 液水, 水蒸気以外の水物質混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real :: tmpq, pow1, pow2, pow3, coe1, coe2, coe3, epsi, tmpv, tmpcoe

  epsi=Rd/Rv
  tmpv=qvT_2_Tv( qv, T )

  if(present(qo))then
     n=size(qo)
     tmpq=ql+qv
     do i=1,n
        tmpq=tmpq+qo(i)
     end do
  else
     tmpq=ql+qv
  end if

  pow1=(Rd+tmpq*Rv)/(Cpd+tmpq*Cpv)
  pow2=(tmpq*Rv)/(Cpd+tmpq*Cpv)
  pow3=-LH( T )*ql/((Cpd+tmpq*Cpv)*T)
  tmpcoe=(p0/pres)**pow1
  coe1=(1.0-ql/(epsi+tmpq))**(pow1-1.0)
  coe2=(tmpq/(tmpq-ql))**pow2
  coe3=exp(pow3)

  thetalv_Emanuel=tmpv*tmpcoe*coe1*coe2*coe3

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetalv_Emanuel', thetalv_Emanuel, 'K' )
  end if

  return
end function

real function thetavqvP_2_T( ptv, qv, P, dl )  ! 仮温位, 水蒸気混合比, 圧力から温度を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: ptv ! 仮温位 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg / kg]
  real, intent(in) :: P   ! 圧力 [Pa]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: kappa, Tv

  kappa=Rd/cpd
  Tv=ptv*exner_func_dry( P )
  thetavqvP_2_T=qvTv_2_T( qv, Tv )

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetavqvP_2_T', thetavqvP_2_T, 'K' )
  end if

  return
end function

real function thetaw_Emanuel( T, qv, pres, eps, dl )
! Emanuel (1994) の (4.7.10) で定義される湿球温位を計算する.
  use Thermo_Const
  implicit none
  real, intent(in) :: T  ! 温度 [K]
  real, intent(in) :: qv  ! 水蒸気混合比 [kg/kg]
  real, intent(in) :: pres  ! 圧力 [Pa]
  integer, intent(in), optional :: eps  ! 収束条件 (デフォルト = 1.e-6)
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, n
  real, parameter :: a=0.81, b=3376.0, c=2.54, theta1=100.0, theta2=600.0
  real :: tmpa, tmpb, tmpc, fa, fb, fc, tmpthetae, thetaa, thetab, thetac
  real :: rast, err_max, tmp_err

  if(present(eps))then
     err_max=eps
  else
     err_max=1.0e-6
  end if

  tmpthetae=thetae_Bolton( T, qv, pres )
  rast=TP_2_qvs( T, pres )

  tmp_err=err_max

  tmpa=rast*(1.0+a*rast)*(b/theta1-c)-log(tmpthetae/theta1)
  tmpb=rast*(1.0+a*rast)*(b/theta2-c)-log(tmpthetae/theta2)
  thetaa=theta1
  thetab=theta2

  do while (err_max<=tmp_err)
     thetac=0.5*(thetab+thetaa)
     tmpc=rast*(1.0+a*rast)*(b/thetac-c)-log(tmpthetae/thetac)

     if(tmpa*tmpc>0.0)then
        tmp_err=abs(tmpa-tmpc)
        tmpa=tmpc
        thetaa=thetac
     else
        tmp_err=abs(tmpb-tmpc)
        tmpb=tmpc
        thetab=thetac
     end if
  end do

  thetaw_Emanuel=thetac

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Function', 'thetaw_Emanuel', thetaw_Emanuel, 'K' )
  end if

  return
end function