!c Description:
!c  AMP 法による飽和蒸気圧, 飽和蒸気圧の 1 回, 2 回微分
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-02-06    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_svap_amp
  interface

     subroutine svap_amp(temp, svap_cef, svap, svap_dt, svap_2dt, type)  
       use prm_type, only: type_svap
       real(8), intent(in)            :: temp      !温度
       real(8), intent(out)           :: svap      !飽和蒸気圧
       real(8), intent(out), optional :: svap_dt   !飽和蒸気圧の温度微分
       real(8), intent(out), optional :: svap_2dt  !飽和蒸気圧の温度の2階微分
       type(type_svap), intent(in)    :: svap_cef  !飽和蒸気圧式の係数 
       integer, intent(out), optional :: type
     end subroutine svap_amp

  end interface
end module interface_svap_amp


subroutine svap_amp(temp, svap_cef, svap, svap_dt, svap_2dt, type)  
  use prm_phys, only: gasr
  use prm_type, only: type_svap
  
  implicit none
  real(8), intent(in)            :: temp         !温度
  real(8), intent(out)           :: svap         !飽和蒸気圧
  real(8), intent(out), optional :: svap_dt      !飽和蒸気圧の温度微分
  real(8), intent(out), optional :: svap_2dt     !飽和蒸気圧の温度の2階微分
  type(type_svap), intent(in)    :: svap_cef     !飽和蒸気圧式の係数  
  integer, intent(out), optional :: type

  real(8)                        :: svap_log     !飽和蒸気圧の log
  real(8)                        :: svap_log_dt  !飽和蒸気圧の log の温度微分
  real(8)                        :: svap_log_2dt !飽和蒸気圧のlogの2階微分


!!!
!!!飽和蒸気圧を求める
!!!
  svap_log = &
       & svap_cef%A / temp &
       & + svap_cef%B &
       & + (svap_cef%C * dlog(temp)) &
       & + (svap_cef%D * temp) &
       & + (svap_cef%E * (temp**2)) &
       & + dlog(1.0D-1)
  
  svap_log_dt = &
       & - (svap_cef%A / temp**2.0D0) &
       & + (svap_cef%C / temp) &
       & + (svap_cef%D) &
       & + (2.0D0 * svap_cef%E * temp) 
  
  svap_log_2dt = &
       & (2.0D0 * svap_cef%A / temp**3.0D0) &
       & - (svap_cef%C / temp**2.0D0) &
       & + (2.0D0 * svap_cef%E)

!!!
!!!出力する変数を準備
!!!温度範囲をチェック. 欠損値は 999.
!!!
  if (temp < svap_cef%low%temp) then 
     svap = 999.0d0
     if (present(svap_dt)) svap_dt = 999.0d0
     if (present(svap_2dt)) svap_2dt = 999.0d0     
     if (present(type)) type = 1
  elseif (temp > svap_cef%high%temp) then 
     svap = 999.0d0
     if (present(svap_dt)) svap_dt = 999.0d0
     if (present(svap_2dt)) svap_2dt = 999.0d0
     if (present(type)) type = 2
  else
     svap = dexp(svap_log)
     if (present(svap_dt)) svap_dt = dexp(svap_log_dt)
     if (present(svap_2dt)) svap_2dt = dexp(svap_log_2dt)
     if (present(type)) type = 0
  end if

  
  
end subroutine svap_amp
