real function CAPE( tmp_p, tmp_z, tmp_qv, tmp_t, z_ref, undeff, opt, dl )
! 対流有効位置エネルギー (Convective Available Potential Energy) を計算する.
  use Thermo_Const
  use Thermo_Function
  use Algebra
  use Statistics
  implicit none
  real, intent(in) :: tmp_p(:)  ! 気圧 [Pa]
  real, intent(in) :: tmp_z(size(tmp_p))  ! 高度 [m]
  real, intent(in) :: tmp_qv(size(tmp_p))  ! 混合比 [kg/kg]
  real, intent(in) :: tmp_t(size(tmp_p))  ! 温度 [K]
  real, intent(in) :: z_ref  ! パーセルを持ち上げる基準高度 [m]
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: opt  ! LNB 計算の際のオプション (z_LNB 参照)
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: nx, i, iz_LFC, iz_LNB, option
  real :: PLFC, TLFC, PLNB, ZLFC, ZLNB
  real :: tmp(size(tmp_p)), t_par(size(tmp_p))
  real :: inv_tmp(size(tmp_p)), inv_p(size(tmp_p))
  real :: undef

  t_par=0.0

  nx=size(tmp_p)

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

  if(present(opt))then
     option=opt
  else
     option=1
  end if

!-- LFC, LNB での高度, 温度, 圧力の計算

  TLFC=T_LFC( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv )
  ZLFC=z_LFC( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv, undeff=undef )
  ZLNB=z_LNB( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv, opt=option, undeff=undef )

  if(ZLFC/=undef.and.ZLNB/=undef)then
     call interpo_search_1d( tmp_z, ZLFC, iz_LFC, int(undef) )
     call interpo_search_1d( tmp_z, ZLNB, iz_LNB, int(undef) )
     if(iz_LFC>=nx)then
        iz_LFC=int(undef)
     end if
     if(iz_LNB>=nx)then
        iz_LNB=int(undef)
     end if
  else
     iz_LFC=int(undef)
     iz_LNB=int(undef)
  end if

  if(iz_LFC/=int(undef).and.iz_LNB/=int(undef))then

  if(tmp_t(iz_LFC)/=undef.and.tmp_qv(iz_LFC)/=undef.and.tmp_p(iz_LFC)/=undef .and.tmp_t(iz_LFC+1)/=undef.and.tmp_qv(iz_LFC+1)/=undef.and.tmp_p(iz_LFC+1)/=undef .and.tmp_t(iz_LNB)/=undef.and.tmp_qv(iz_LNB)/=undef.and.tmp_p(iz_LNB)/=undef .and.tmp_t(iz_LNB+1)/=undef.and.tmp_qv(iz_LNB+1)/=undef.and.tmp_p(iz_LNB+1)/=undef)then
     call interpolation_1d( (/tmp_z(iz_LFC), tmp_z(iz_LFC+1)/), (/tmp_p(iz_LFC), tmp_p(iz_LFC+1)/), ZLFC, PLFC )
     call interpolation_1d( (/tmp_z(iz_LNB), tmp_z(iz_LNB+1)/), (/tmp_p(iz_LNB), tmp_p(iz_LNB+1)/), ZLNB, PLNB )

     do i=1,nx
        if(iz_LFC<=i.and.i<=iz_LNB+1)then

           if(undef==tmp_t(i))then
              tmp(i)=undef
           else
              t_par(i)=moist_laps_temp( PLFC, TLFC, tmp_p(i) )
              tmp(i)=(t_par(i)-tmp_t(i))/tmp_p(i)
           end if
        else
           tmp(i)=undef
        end if
     end do

!-- 積分のため, 順序の入れ替え
     do i=1,nx
        inv_tmp(i)=tmp(nx-i+1)
        inv_p(i)=tmp_p(nx-i+1)
     end do

     call rectangle_int( inv_p, inv_tmp, PLNB, PLFC, CAPE, undef )

     CAPE=CAPE*Rd

  else

     CAPE=undef

  end if

  else

     CAPE=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'CAPE', CAPE, 'J kg-1' )
  end if

  return
end function

real function CIN( tmp_p, tmp_z, tmp_qv, tmp_t, z_ref, undeff, dl )
! 対流抑制エネルギー (Convective INhibitation energy) を計算する.
  use Thermo_Const
  use Algebra
  use Thermo_Function
  use Statistics
  implicit none
  real, intent(in) :: tmp_p(:)  ! 気圧 [Pa]
  real, intent(in) :: tmp_z(size(tmp_p))  ! 高度 [m]
  real, intent(in) :: tmp_qv(size(tmp_p))  ! 混合比 [kg/kg]
  real, intent(in) :: tmp_t(size(tmp_p))  ! 温度 [K]
  real, intent(in) :: z_ref  ! パーセルを持ち上げる基準高度 [m]
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: nx, i, iz_LFC, iz_ref, iz_LCL
  real :: PLFC, TLFC, ZLFC, PLCL, ZLCL, TLCL, p_ref, T_ref, qv_ref
  real :: tmp(size(tmp_p)), t_par(size(tmp_p))
  real :: inv_tmp(size(tmp_p)), inv_p(size(tmp_p))
  real :: undef

  t_par=0.0

  nx=size(tmp_p)

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

!-- z_ref を基準とした, LCL, LFC での温度, 高度, 圧力の計算

!-- LFC 高度の変数計算
  TLFC=T_LFC( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv )
  ZLFC=z_LFC( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv )

  call interpo_search_1d( tmp_z, ZLFC, iz_LFC, int(undef) )
  call interpo_search_1d( tmp_z, z_ref, iz_ref, int(undef) )

  if(iz_LFC>=nx)then
     iz_LFC=int(undef)
  end if
  if(iz_ref>=nx)then
     iz_ref=int(undef)
  end if

  if(iz_LFC/=int(undef).and.iz_ref/=int(undef))then
  if(tmp_t(iz_LFC)/=undef.and.tmp_qv(iz_LFC)/=undef.and.tmp_p(iz_LFC)/=undef .and.tmp_t(iz_LFC+1)/=undef.and.tmp_qv(iz_LFC+1)/=undef.and.tmp_p(iz_LFC+1)/=undef .and.tmp_t(iz_ref)/=undef.and.tmp_qv(iz_ref)/=undef.and.tmp_p(iz_ref)/=undef .and.tmp_t(iz_ref+1)/=undef.and.tmp_qv(iz_ref+1)/=undef.and.tmp_p(iz_ref+1)/=undef)then

     call interpolation_1d( (/tmp_z(iz_LFC), tmp_z(iz_LFC+1)/), (/tmp_p(iz_LFC), tmp_p(iz_LFC+1)/), ZLFC, PLFC )
!-- reference 高度の変数計算
     call interpolation_1d( (/tmp_z(iz_ref), tmp_z(iz_ref+1)/), (/tmp_p(iz_ref), tmp_p(iz_ref+1)/), z_ref, p_ref )
     call interpolation_1d( (/tmp_z(iz_ref), tmp_z(iz_ref+1)/), (/tmp_t(iz_ref), tmp_t(iz_ref+1)/), z_ref, T_ref )
     call interpolation_1d( (/tmp_z(iz_ref), tmp_z(iz_ref+1)/), (/tmp_qv(iz_ref), tmp_qv(iz_ref+1)/), z_ref, qv_ref )
!-- LCL 高度の変数計算
     TLCL=TqvP_2_TLCL( t_ref, qv_ref, p_ref )
     ZLCL=z_LCL( z_ref, tmp_z, tmp_t, tmp_p, tmp_qv )

     call interpo_search_1d( tmp_z, ZLCL, iz_LCL, int(undef) )

     if(iz_LCL>=nx)then
        iz_LCL=int(undef)
     end if

     if(iz_LCL/=int(undef))then
     if(tmp_t(iz_LCL)/=undef.and.tmp_qv(iz_LCL)/=undef.and.tmp_p(iz_LCL)/=undef .and.tmp_t(iz_LCL+1)/=undef.and.tmp_qv(iz_LCL+1)/=undef.and.tmp_p(iz_LCL+1)/=undef)then

        call interpolation_1d( (/tmp_z(iz_LCL), tmp_z(iz_LCL+1)/), (/tmp_p(iz_LCL), tmp_p(iz_LCL+1)/), ZLCL, PLCL )

! 周囲の温度は T で得られているので, まず LCL を計算し, その高度
! までを乾燥断熱で, それ以降を湿潤断熱でパーセル温度を計算.

        do i=1,nx
           if(i>=iz_ref.and.i<=iz_LCL+1)then
              t_par(i)=T_ref-(g/Cpd)*(tmp_z(i)-z_ref)
           else
              if(i>iz_LCL+1.and.i<=iz_LFC+1)then
                 t_par(i)=moist_laps_temp( PLCL, TLCL, tmp_p(i) )
              else
                 t_par(i)=undef
              end if
           end if

           if(undef==tmp_t(i).or.t_par(i)==undef)then
              tmp(i)=undef
           else
              tmp(i)=(t_par(i)-tmp_t(i))/tmp_p(i)
           end if
        end do

!do i=1,nx
!   if(tmp_p(i)<p_ref.and.tmp_p(i)>PLFC)then
!      t_par(i)=t_bot-(g/CPD)*(height(i+1)-height(i))
!      if(PLCL>PLFC)then   ! 強い逆転層のある場合
!         if(p(i)<=PLCL.and.p(i)>=PLFC)then
!            call moist_laps_calc( PLCL, TLCL, p(i), t_par(i) )
!         end if
!      end if
!   end if
!end do

!-- 積分のため, 順序の入れ替え
        do i=1,nx
           inv_tmp(i)=tmp(nx-i+1)
           inv_p(i)=tmp_p(nx-i+1)
        end do

        call rectangle_int( inv_p, inv_tmp, PLFC, p_ref, CIN, undef )

        CIN=CIN*Rd

     else

        CIN=undef

     end if

     else

        CIN=undef

     end if

  else

     CIN=undef

  end if

  else

     CIN=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'CIN', CIN, 'J kg-1' )
  end if

  return
end function

real function Louis( z, z0m, richard, dl )
! Louis(1980) で提案されている大気の不安定度を考慮したバルク係数の補正係数を計算する関数
  use Thermo_Const
  implicit none
  real, intent(in) :: z  ! cm を求める高度 [m]
  real, intent(in) :: z0m  ! モデルで計算される粗度高度 [m]
  real, intent(in) :: richard  ! バルクリチャードソン数
  integer, intent(in), optional :: dl  ! デバッグレベル
  real, parameter :: b=5.0, c=5.0
  real :: cm_tmp, zratio

  cm_tmp=(kalm/(log(z)-log(z0m)))**2
  zratio=z/z0m

  if(richard<0.0)then
     Louis=1.0-((2.0*b*richard)/(1.0+3.0*b*c*cm_tmp*sqrt(-richard*zratio)))
  else
     Louis=1.0/(1.0+2.0*b*richard*sqrt(1.0+c*richard))
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'Louis', Louis, '1' )
  end if

  return
end function

real function Rich( za, pta, ptg, va, qva, qvs, dl )
! バルクリチャードソン数を計算する関数
  use Phys_Const
  use Thermo_Function
  implicit none
  real, intent(in) :: za  ! リチャードソン数を計算する高度 [m]
  real, intent(in) :: pta  ! za での仮温位 [K]
  real, intent(in) :: ptg  ! 地表面での温位 [K]
  real, intent(in) :: va  ! 高度 za での水平風速の絶対値 [m/s]
  real, intent(in) :: qva  ! za での混合比 [kg/kg]
  real, intent(in) :: qvs  ! 地表面での飽和混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: ptvg, ptva, dpt

  ptvg=ptg*((1.0+eps_rdrv*qvs)/(1.0+qvs))
  ptva=pta*((1.0+eps_rdrv*qva)/(1.0+qva))
  dpt=ptva-ptvg
  Rich=(g*za*dpt)/(ptva*(va**2))

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'Rich', Rich, '1' )
  end if

  return
end function

real function T_LFC( z_ref, z, temp, pres, qv, undeff, dl )  ! LFC 高度での温度を計算する.
! ある基準高度での相当温位を上昇させたとき, 周囲の飽和相当温位と一致する高度が
! LFC である. 
! ここでは, 周囲の飽和相当温位と交差する 2 点を求め, そこから高度方向に
! 内挿することで高度を決定する.
  use Statistics
  use Thermo_Function
  implicit none
  real, intent(in) :: z_ref  ! 基準高度 [m]
  real, intent(in) :: z(:)  ! 高度座標 [m]
  real, intent(in) :: temp(size(z))  ! 温度 [K]
  real, intent(in) :: pres(size(z))  ! 圧力 [Pa]
  real, intent(in) :: qv(size(z))  ! 水蒸気混合比 [kg/kg]
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, nz, iz, iept
  real :: sept(size(z))
  real :: eptiz, eptiz1, ept_ref, z_sol
  real :: undef

  nz=size(z)
  iept=0

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

!-- 鉛直方向の飽和相当温位を計算.
  do i=1,nz
     sept(i)=thetaes_Bolton( temp(i), pres(i) )
  end do

!-- z_ref での相当温位を計算する.
  call interpo_search_1d( z, z_ref, iz, int(undef) )
  if(iz>=nz)then
     iz=int(undef)
  end if

!-- iz と iz+1 での相当温位を計算.
  if(iz/=int(undef))then
  if(temp(iz)/=undef.and.qv(iz)/=undef.and.pres(iz)/=undef .and.temp(iz+1)/=undef.and.qv(iz+1)/=undef.and.pres(iz+1)/=undef)then
     eptiz=thetae_Bolton( temp(iz), qv(iz), pres(iz) )
     eptiz1=thetae_Bolton( temp(iz+1), qv(iz+1), pres(iz+1) )
!-- これらの間から相当温位を内挿
     call interpolation_1d( (/z(iz), z(iz+1)/), (/eptiz, eptiz1/), z_ref, ept_ref )
!-- z_ref より上の高度で ept_ref と飽和相当温位の関係を計算.
     do i=iz+1,nz-1
        if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
           iept=i
           exit
        end if
     end do
!-- パーセルの相当温位が環境場の飽和相当温位と逆転する直前の高度が iept である.
!-- ここで, iept と iept+1 の飽和相当温位から高度を内挿する.
     if(iept>0.and.iept<nz)then
        call interpolation_1d( (/sept(iept), sept(iept+1)/), (/temp(iept), temp(iept+1)/), ept_ref, z_sol )

        T_LFC=z_sol
     else
        T_LFC=undef
     end if

  else

     T_LFC=undef

  end if

  else

     T_LFC=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'T_LFC', T_LFC, 'K' )
  end if

  return
end function

real function T_LNB( z_ref, z, temp, pres, qv, opt, undeff, dl )  ! LNB 高度での温度を計算する.
! LFC 高度以上で再び環境場の飽和相当温位と交わる高度が LNB である.
! ただし, 実際の測定データでは, 測器の微小変動により, LFC の直上で LNB に
! 到達する場合がある.
! そこで, opt として, LFC を求めてから LNB を計算する方法と
! 上空から下向きに下ろして最初に交差する高度を LNB と定義する手法の
! 2 パターン用意することにする.
! opt = 1, LFC から求める. opt = 2 直上から計算.
! デフォルトでは opt = 1 を計算.
  use Statistics
  use Thermo_Function
  implicit none
  real, intent(in) :: z_ref  ! 基準高度 [m]
  real, intent(in) :: z(:)  ! 高度座標 [m]
  real, intent(in) :: temp(size(z))  ! 温度 [K]
  real, intent(in) :: pres(size(z))  ! 圧力 [Pa]
  real, intent(in) :: qv(size(z))  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: opt  ! 計算手法のオプション
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, nz, iz, iept, counter
  real :: sept(size(z))
  real :: eptiz, eptiz1, ept_ref, z_sol
  real :: undef

  nz=size(z)
  iept=0

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

!-- 鉛直方向の飽和相当温位を計算.
  do i=1,nz
     sept(i)=thetaes_Bolton( temp(i), pres(i) )
  end do

!-- z_ref での相当温位を計算する.
  call interpo_search_1d( z, z_ref, iz, int(undef) )
  if(iz>=nz)then
     iz=int(undef)
  end if

!-- iz と iz+1 での相当温位を計算.
  if(iz/=int(undef))then
  if(temp(iz)/=undef.and.qv(iz)/=undef.and.pres(iz)/=undef .and.temp(iz+1)/=undef.and.qv(iz+1)/=undef.and.pres(iz+1)/=undef)then
     eptiz=thetae_Bolton( temp(iz), qv(iz), pres(iz) )
     eptiz1=thetae_Bolton( temp(iz+1), qv(iz+1), pres(iz+1) )
!-- これらの間から相当温位を内挿
     call interpolation_1d( (/z(iz), z(iz+1)/), (/eptiz, eptiz1/), z_ref, ept_ref )

     counter=0

     if(present(opt))then
        if(opt==2)then
           do i=nz,iz+1,-1  ! 上から下に下ろす.
              if((sept(i)-ept_ref)*(sept(i-1)-ept_ref)<0.0)then
                 iept=i-1  ! 上から下に下ろしているので, 1 要素下のデータが iept.
                 exit
              end if
           end do
        else
!-- z_ref より上の高度で ept_ref と飽和相当温位の関係を計算.
           do i=iz+1,nz-1
              if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
                 counter=counter+1
                 if(counter==2)then
                    iept=i
                    exit
                 end if
              end if
           end do
        end if
     else
        do i=iz+1,nz-1
           if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
              counter=counter+1
              if(counter==2)then
                 iept=i
                 exit
              end if
           end if
        end do
     end if
!-- パーセルの相当温位が環境場の飽和相当温位と逆転する直前の高度が iept である.
!-- ここで, iept と iept+1 の飽和相当温位から高度を内挿する.
     if(iept>0.and.iept<nz)then
        call interpolation_1d( (/sept(iept), sept(iept+1)/), (/temp(iept), temp(iept+1)/), ept_ref, z_sol )

        T_LNB=z_sol
     else
        T_LNB=undef
     end if

  else

     T_LNB=undef

  end if

  else

     T_LNB=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'T_LNB', T_LNB, 'K' )
  end if

  return
end function

real function cm( z, z0m, richard, dl )
! 運動量に関するバルク係数を計算する関数
  use Thermo_Const
  implicit none
  real, intent(in) :: z  ! cm を求める高度 [m]
  real, intent(in) :: z0m  ! モデルで計算される粗度高度 [m]
  real, intent(in), optional :: richard  ! Louis (1980) のスキームで計算する場合のバルクリチャードソン数
  integer, intent(in), optional :: dl  ! デバッグレベル

  if(present(richard))then
     cm=(kalm/(log(z)-log(z0m)))**2
     cm=cm*Louis( z, z0m, richard )
  else
     cm=(kalm/(log(z)-log(z0m)))**2
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'cm', cm, '1' )
  end if

  return
end function

real function cmdva_2_ust( cmd, va, dl )
! バルク係数, 速度から摩擦速度 u_* を計算する関数
  implicit none
  real, intent(in) :: cmd  ! 高度 za でのバルク係数
  real, intent(in) :: va  ! 高度 za での水平風の絶対値 [m/s]
  integer, intent(in), optional :: dl  ! デバッグレベル

  cmdva_2_ust=va*sqrt(cmd)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'cmdva_2_ust', cmdva_2_ust, 'm s-1' )
  end if

end function

real function precip_water( p, qv, undef, dl )  ! 可降水量を計算する. 単位は [kg/m^2]
! 積分範囲は p の要素数の上下端で自動的に指定.
! 気象がわかる数と式より, 高度座標で積分するのではなく,
! 静力学平衡から気圧座標へ置き直して積分.
  use Algebra
  use Phys_Const
  implicit none
  real, intent(in) :: p(:)  ! 圧力 [Pa]
  real, intent(in) :: qv(size(p))  ! 水蒸気混合比 [kg/kg]
  real, intent(in), optional :: undef  ! undef
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: nx, i
  real, dimension(size(p)) :: tmp_p, tmp_qv
  real :: precip

  nx=size(p)

!-- 積分のため, 順序の入れ替え
  do i=1,nx
     tmp_qv(i)=qv(nx-i+1)
     tmp_p(i)=p(nx-i+1)
  end do

  if(present(undef))then
     call rectangle_int( tmp_p, tmp_qv, tmp_p(1), tmp_p(nx), precip, undef )
  else
     call rectangle_int( tmp_p, tmp_qv, tmp_p(1), tmp_p(nx), precip )
  end if

  precip_water=precip/g

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'precip_water', precip_water, 'mm' )
  end if

  return
end function

real function qrsg_2_dbz( rho, qr, qs, qg, dl )
! 雨水混合比から擬似的なレーダー反射強度に変換する関数.
! レーダー反射強度 dbz はレーダー反射因子 Z を用いて,
! $$ dbz = 10 \log{Z} $$
! で定義される. さらに, 反射因子 Z は
! $$Z=\sum_{D^6} =int^{\infty}_{0}{N(D)D^6dD} $$
! という式で定義される. ここで, D は水物質の直径, N は数密度である.
! ここでは気象レーダーを仮定しているので, 降水粒子のカテゴリは
! 雨水, 雪, 霰 のみとする.
! warm rain の場合でも, 雪と霰の混合比をゼロで引数に入れることで計算が可能.
! 粒径分布は指数関数分布を仮定して計算する. その場合, 各カテゴリ粒子の粒径分布は
! 粒子の直径を D_x として, 
! $$n_x(D_x)=n_{x0}\exp{(-\lambda _xD_x)} $$
! という式で表現される. ここで, $\lambda _x$は粒径分布の傾きである.
! 変換公式としては, Murakami (1990) の (54) 式を用いて計算を行う.
  use Math_Const
  use Cloud_Const
  use Special_Function
  implicit none
  real, intent(in) :: rho  ! 大気の密度 [kg/m^3]
  real, intent(in) :: qr  ! 雨水の混合比 [kg/kg]
  real, intent(in), optional :: qs  ! 雪の混合比 [kg/kg]
  real, intent(in), optional :: qg  ! 霰の混合比 [kg/kg]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: coef
  real :: z_tmp

  z_tmp=0.0
  coef=kaijo(6.0)*(((epi-1.0)/(epi+2.0))/((epw-1.0)/(epw+2.0))**2)

  if(present(qs))then
     z_tmp=z_tmp+coef*((rhos/rhow)**2)*ns0*((rho*qs/(pi*rhos*ns0))**1.75)!*1.0e18
  end if

  if(present(qg))then
     z_tmp=z_tmp+coef*((rhog/rhow)**2)*ng0*((rho*qg/(pi*rhog*ng0))**1.75)!*1.0e18
  end if

  z_tmp=z_tmp+nr0*((rho*qr/(pi*rhow*nr0))**1.75)!*1.0e18
  z_tmp=coef*z_tmp

  if(z_tmp>0.0)then
     qrsg_2_dbz=180.0+10.0*log10(z_tmp)
  else
     qrsg_2_dbz=0.0
  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'qrsg_2_dbz', qrsg_2_dbz, 'dBZ' )
  end if

  return
end function

real function taurho_2_ust( tau, rho, dl )
! 風応力のデカルト水平 2 成分とその高度での密度から摩擦速度 u_* を計算する関数
  implicit none
  real, intent(in) :: tau(2)  ! 基準高度での風応力のデカルト水平 2 成分 [N/m]
  real, intent(in) :: rho  ! 基準高度での密度 [kg/m^3]
  integer, intent(in), optional :: dl  ! デバッグレベル
  real :: taua

  taua=sqrt(tau(1)**2+tau(2)**2)
  taurho_2_ust=sqrt(taua/rho)

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'taurho_2_ust', taurho_2_ust, 'm s-1' )
  end if

end function

real function z_LCL( z_ref, z, temp, pres, qv, undeff, dl )  ! LCL 高度を計算する.
! パーセルが z_ref から断熱過程で上昇し, この温度に達したときの高度が
! LCL 高度となるので, z_ref から乾燥断熱減率で以下のように計算する.
! LCL 高度を z_LCL とし, パーセルの基準高度を z_ref とすると,
! LCL に達するまでは乾燥断熱減率 \Gamma _d で変化するので,
! $\Gamma _d=\frac{g}{C_p} =-\frac{T_LCL-T_ref}{z_LCL-z_ref} $
! という式が成り立つ. ここで, z_LCL について解くと,
! $z_LCL=z_ref+\frac{C_p}{g} (T_ref-T_LCL)$
! となる.
  use Thermo_Const
  use Thermo_Function
  use Phys_Const
  use Statistics
  implicit none
  real, intent(in) :: z_ref  ! 基準高度 [m]
  real, intent(in) :: z(:)  ! 高度座標 [m]
  real, intent(in) :: temp(size(z))  ! 温度 [K]
  real, intent(in) :: pres(size(z))  ! 圧力 [Pa]
  real, intent(in) :: qv(size(z))  ! 水蒸気混合比 [kg/kg]
  real ,intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: iz_ref
  real :: TLCL
  real :: T_ref, P_ref, qv_ref
  real :: undef

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

  call interpo_search_1d( z, z_ref, iz_ref, int(undef) )

  if(iz_ref>=size(z))then
     iz_ref=int(undef)
  end if

  if(iz_ref/=int(undef))then
  if(temp(iz_ref)/=undef.and.qv(iz_ref)/=undef.and.pres(iz_ref)/=undef .and.temp(iz_ref+1)/=undef.and.qv(iz_ref+1)/=undef.and.pres(iz_ref+1)/=undef)then
     call interpolation_1d( (/z(iz_ref), z(iz_ref+1)/), (/temp(iz_ref), temp(iz_ref+1)/), z_ref, T_ref )
     call interpolation_1d( (/z(iz_ref), z(iz_ref+1)/), (/pres(iz_ref), pres(iz_ref+1)/), z_ref, P_ref )
     call interpolation_1d( (/z(iz_ref), z(iz_ref+1)/), (/qv(iz_ref), qv(iz_ref+1)/), z_ref, qv_ref )

     TLCL=TqvP_2_TLCL( T_ref, qv_ref, P_ref )  ! z_ref のパーセルの LCL での温度.
     z_LCL=z_ref+(Cpd/g)*(T_ref-TLCL)

  else

     z_LCL=undef

  end if

  else

     z_LCL=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'z_LCL', z_LCL, 'm' )
  end if

  return
end function

real function z_LFC( z_ref, z, temp, pres, qv, undeff, dl )  ! LFC 高度を計算する.
! ある基準高度での相当温位を上昇させたとき, 周囲の飽和相当温位と一致する高度が
! LFC である. 
! ここでは, 周囲の飽和相当温位と交差する 2 点を求め, そこから高度方向に
! 内挿することで高度を決定する.
  use Statistics
  use Thermo_Function
  implicit none
  real, intent(in) :: z_ref  ! 基準高度 [m]
  real, intent(in) :: z(:)  ! 高度座標 [m]
  real, intent(in) :: temp(size(z))  ! 温度 [K]
  real, intent(in) :: pres(size(z))  ! 圧力 [Pa]
  real, intent(in) :: qv(size(z))  ! 水蒸気混合比 [kg/kg]
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, nz, iz, iept
  real :: sept(size(z))
  real :: eptiz, eptiz1, ept_ref, z_sol
  real :: undef

  nz=size(z)
  iept=0

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

!-- 鉛直方向の飽和相当温位を計算.
  do i=1,nz
     sept(i)=thetaes_Bolton( temp(i), pres(i) )
  end do

!-- z_ref での相当温位を計算する.
  call interpo_search_1d( z, z_ref, iz, int(undef) )

  if(iz>=nz)then
     iz=int(undef)
  end if

  if(iz/=int(undef))then
  if(temp(iz)/=undef.and.qv(iz)/=undef.and.pres(iz)/=undef .and.temp(iz+1)/=undef.and.qv(iz+1)/=undef.and.pres(iz+1)/=undef)then
!-- iz と iz+1 での相当温位を計算.
     eptiz=thetae_Bolton( temp(iz), qv(iz), pres(iz) )
     eptiz1=thetae_Bolton( temp(iz+1), qv(iz+1), pres(iz+1) )
!-- これらの間から相当温位を内挿
     call interpolation_1d( (/z(iz), z(iz+1)/), (/eptiz, eptiz1/), z_ref, ept_ref )
!-- z_ref より上の高度で ept_ref と飽和相当温位の関係を計算.
     do i=iz+1,nz-1
        if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
           iept=i
           exit
        end if
     end do
!-- パーセルの相当温位が環境場の飽和相当温位と逆転する直前の高度が iept である.
!-- ここで, iept と iept+1 の飽和相当温位から高度を内挿する.
     if(iept>0.and.iept<nz)then
        call interpolation_1d( (/sept(iept), sept(iept+1)/), (/z(iept), z(iept+1)/), ept_ref, z_sol )

        z_LFC=z_sol
     else
        z_LFC=undef
     end if

  else

     z_LFC=undef

  end if

  else

     z_LFC=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'z_LFC', z_LFC, 'm' )
  end if

  return
end function

real function z_LNB( z_ref, z, temp, pres, qv, opt, undeff, dl )  ! LNB 高度を計算する.
! LFC 高度以上で再び環境場の飽和相当温位と交わる高度が LNB である.
! ただし, 実際の測定データでは, 測器の微小変動により, LFC の直上で LNB に
! 到達する場合がある.
! そこで, opt として, LFC を求めてから LNB を計算する方法と
! 上空から下向きに下ろして最初に交差する高度を LNB と定義する手法の
! 2 パターン用意することにする.
! opt = 1, LFC から求める. opt = 2 直上から計算.
! デフォルトでは opt = 1 を計算.
  use Statistics
  use Thermo_Function
  implicit none
  real, intent(in) :: z_ref  ! 基準高度 [m]
  real, intent(in) :: z(:)  ! 高度座標 [m]
  real, intent(in) :: temp(size(z))  ! 温度 [K]
  real, intent(in) :: pres(size(z))  ! 圧力 [Pa]
  real, intent(in) :: qv(size(z))  ! 水蒸気混合比 [kg/kg]
  integer, intent(in), optional :: opt  ! 計算手法のオプション
  real, intent(in), optional :: undeff
  integer, intent(in), optional :: dl  ! デバッグレベル
  integer :: i, nz, iz, iept, counter
  real :: sept(size(z))
  real :: eptiz, eptiz1, ept_ref, z_sol
  real :: undef

  nz=size(z)
  iept=0

  if(present(undeff))then
     undef=undeff
  else
     undef=-999.0
  end if

!-- 鉛直方向の飽和相当温位を計算.
  do i=1,nz
     sept(i)=thetaes_Bolton( temp(i), pres(i) )
  end do

!-- z_ref での相当温位を計算する.
  call interpo_search_1d( z, z_ref, iz, int(undef) )

  if(iz>=nz)then
     iz=int(undef)
  end if

  if(iz/=int(undef))then
  if(temp(iz)/=undef.and.qv(iz)/=undef.and.pres(iz)/=undef .and.temp(iz+1)/=undef.and.qv(iz+1)/=undef.and.pres(iz+1)/=undef)then
!-- iz と iz+1 での相当温位を計算.
     eptiz=thetae_Bolton( temp(iz), qv(iz), pres(iz) )
     eptiz1=thetae_Bolton( temp(iz+1), qv(iz+1), pres(iz+1) )
!-- これらの間から相当温位を内挿
     call interpolation_1d( (/z(iz), z(iz+1)/), (/eptiz, eptiz1/), z_ref, ept_ref )

     counter=0

     if(present(opt))then
        if(opt==2)then
           do i=nz,iz+1,-1  ! 上から下に下ろす.
              if((sept(i)-ept_ref)*(sept(i-1)-ept_ref)<0.0)then
                 iept=i-1  ! 上から下に下ろしているので, 1 要素下のデータが iept.
                 exit
              end if
           end do
        else
!-- z_ref より上の高度で ept_ref と飽和相当温位の関係を計算.
           do i=iz+1,nz-1
              if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
                 counter=counter+1
                 if(counter==2)then
                    iept=i
                    exit
                 end if
              end if
           end do
        end if
     else
        do i=iz+1,nz-1
           if((sept(i)-ept_ref)*(sept(i+1)-ept_ref)<0.0)then
              counter=counter+1
              if(counter==2)then
                 iept=i
                 exit
              end if
           end if
        end do
     end if
!-- パーセルの相当温位が環境場の飽和相当温位と逆転する直前の高度が iept である.
!-- ここで, iept と iept+1 の飽和相当温位から高度を内挿する.
     if(iept>0.and.iept<nz)then
        call interpolation_1d( (/sept(iept), sept(iept+1)/), (/z(iept), z(iept+1)/), ept_ref, z_sol )

        z_LNB=z_sol
     else
        z_LNB=undef
     end if

  else

     z_LNB=undef

  end if

  else

     z_LNB=undef

  end if

  if(present(dl))then
     call debug_flag_r( dl, 'Thermo_Advanced_Function', 'z_LNB', z_LNB, 'm' )
  end if

  return
end function