!c Description:
!c   飽和蒸気圧の経験式の係数を取得
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-02-05    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_def_svap_cef
  interface

     subroutine def_svap_cef(property, propertydir)
       use prm_type, only: type_property
       character(*), intent(in)           :: propertydir
       type(type_property), intent(inout) :: property(:)
     end subroutine def_svap_cef
     
  end interface
end module interface_def_svap_cef



subroutine def_svap_cef(property, propertydir)
  use prm_phys, only: pres_ref, gasr
  use prm_type, only: type_property
  use interface_netcdf_open 
  use interface_netcdf_get
  use interface_netcdf_close
  use interface_svap_amp
  use interface_svap_antoine
  use interface_gibbs_energy
  use interface_entropy
  use interface_name2num

  implicit none
  character(*), intent(in)           :: propertydir
  type(type_property), intent(inout) :: property(:)
  
  real(8)                 :: temp
  real(8)                 :: svap
  real(8)                 :: dt
  real(8)                 :: entr0_dt
  real(8)                 :: entr0_tmp
  real(8)                 :: entr0_tmp_pre
  real(8)                 :: mu0_tmp
  real(8), allocatable    :: rvar(:) 
  real(8), allocatable    :: mol(:) 
  real(8), allocatable    :: entr0(:) 
  real(8), allocatable    :: entr0_pre(:) 
  real(8), allocatable    :: mu0(:) 
  real(8), parameter      :: pres_low = 1.0d1   !1.0d2
  real(8), parameter      :: pres_high = 1.0d20   !1.0d5, 1.0d6(2005-05-07)
  real(8), parameter      :: temp_low = 5.0d1
  real(8), parameter      :: temp_high = 1.0d3
  integer                 :: ncid
  integer                 :: N, s, i
  integer                 :: num1, num2
  integer                 :: num1_sol, num1_liq, num1_cnd
  integer                 :: num2_sol, num2_liq, num2_cnd
  integer, parameter      :: num = 10000
  character(80)           :: name
  character(80)           :: file
  

  !初期化
  N = size(property, 1)
  allocate(mol(N))
  mol = 1.0d0   !ダミー
  property%svap_cef%low%temp = temp_low * 1.0d-1       
  property%svap_cef%low%A = 0.0d0
  property%svap_cef%low%B = 0.0d0
  property%svap_cef%low%C = 0.0d0
  property%svap_cef%high%temp = temp_high * 1.0d1
  property%svap_cef%high%A = 0.0d0
  property%svap_cef%high%B = 0.0d0
  property%svap_cef%high%C = 0.0d0
  property%svap_cef%A = 0.0d0
  property%svap_cef%B = 0.0d0
  property%svap_cef%C = 0.0d0
  property%svap_cef%D = 0.0d0
  property%svap_cef%E = 0.0d0


  calc_svap: do s = 1, N
     !計算方法が 3, 4 以外は終了
     if (property(s)%method == 1 .OR. property(s)%method == 2) cycle calc_svap
     
     !ファイル名準備
     file = trim(PropertyDir) // trim(property(s)%name) // '.nc'
     
     !ファイルオープン
     call netcdf_open(trim(file), ncid)
     
     !飽和蒸気圧
     if (property(s)%method == 3) then 
        name = "amp"
     elseif (property(s)%method == 4) then 
        name = "antoine"
     end if
     call netcdf_get(ncid, name, rvar)

     !ファイルクローズ
     call netcdf_close(ncid)

     !値の代入
     property(s)%svap_cef%A = rvar(1)
     property(s)%svap_cef%B = rvar(2)
     property(s)%svap_cef%C = rvar(3)
     property(s)%svap_cef%D = rvar(4)
     property(s)%svap_cef%E = rvar(5)

!!!
!!! 飽和蒸気圧経験式の適用限界(下限)を計算
!!!
     cri_low: do i = 1, num
        dt = (temp_high - temp_low) / num 
        temp = temp_high - (i-1) * dt
        
        if (property(s)%method == 3) then 
           call svap_amp(temp, property(s)%svap_cef, svap)
        elseif(property(s)%method == 4) then 
           call svap_antoine(temp, property(s)%svap_cef, svap)
        else
           exit cri_low
        end if
        
        if (svap <= pres_low) then 
           call entropy(temp, pres_ref, mol, property, entr0=entr0)
           call gibbs_energy(temp, pres_ref, mol, property, mu0=mu0)
           
           property(s)%svap_cef%low%temp = temp           
           property(s)%svap_cef%low%A = entr0(s)
           property(s)%svap_cef%low%B = mu0(s) + (entr0(s) * temp) / 2.0d0 
           
           deallocate(entr0, mu0)
           exit cri_low           
        end if
     end do cri_low
     
!!!
!!! 飽和蒸気圧経験式の適用限界(上限)を計算
!!!
     cri_high: do i = 1, num
        dt = (temp_high - temp_low) / num 
        temp = temp_high - (i-1) * dt

        !飽和蒸気圧を計算
        if (property(s)%method == 3) then 
           call svap_amp(temp, property(s)%svap_cef, svap)
        elseif(property(s)%method == 4) then 
           call svap_antoine(temp, property(s)%svap_cef, svap)
        else
           exit cri_high
        end if

        !飽和蒸気圧が小さい場合の処理
        if (svap <= pres_high) then 
           call entropy(temp, pres_ref, mol, property, entr0=entr0)
           call entropy(temp - dt, pres_ref, mol, property, entr0=entr0_pre)
           call gibbs_energy(temp, pres_ref, mol, property, mu0=mu0)
           
           entr0_dt = (entr0(s) - entr0_pre(s))/ dt
           property(s)%svap_cef%high%temp = temp                
           property(s)%svap_cef%high%A = entr0_dt
           property(s)%svap_cef%high%B = entr0(s) - entr0_dt * temp
           property(s)%svap_cef%high%C = mu0(s) + entr0(s) * temp &
                & - (entr0_dt * (temp ** 2.0d0)) / 2.0d0 

           deallocate(entr0, entr0_pre, mu0)
           exit cri_high
        end if
     end do cri_high



!!!
!!! NH4SH の扱い
!!!
     !計算番号 5 に対する処置
     if(property(s)%method == 5) then 
        !NH3-g の値を求める
        name = 'NH3-g'
        call name2num(name, property%name, num1)
        name = 'NH3-l'
        call name2num(name, property%name, num1_liq)
        name = 'NH3-s'
        call name2num(name, property%name, num1_sol)
        if (num1_liq /= 0 .AND. num1_liq == 0) then 
           num1_cnd = num1_liq
        elseif (num1_liq == 0 .AND. num1_sol == 0) then 
           write(*,*) "def_svap_cef.f90; No NH3-l, NH3-s"
        else
           num1_cnd = num1_sol
        end if

        !H2S-g の値を求める
        name = 'H2S-g'
        call name2num(name, property%name, num2)
        name = 'H2S-l'
        call name2num(name, property%name, num2_liq)
        name = 'H2S-s'
        call name2num(name, property%name, num2_sol)
        if (num2_liq /= 0 .AND. num2_liq == 0) then 
           num2_cnd = num2_liq
        elseif (num2_liq == 0 .AND. num2_sol == 0) then 
           write(*,*) "def_svap_cef.f90; No H2S-l, H2S-s"
        else
           num2_cnd = num2_sol
        end if


        !下限を求める
        if (property(num1_cnd)%svap_cef%low%temp &
             & <= property(num2_cnd)%svap_cef%low%temp) then 
           temp = property(num2_cnd)%svap_cef%low%temp
        else
           temp = property(num1_cnd)%svap_cef%low%temp
        end if
        
        call entropy(temp, pres_ref, mol, property, entr0=entr0)
        entr0_tmp = entr0(num1) + entr0(num2) &
             & - gasr * 61.781d0  + gasr * dlog(1.0d12) 
        call gibbs_energy(temp, pres_ref, mol, property, mu0=mu0)
        mu0_tmp  = mu0(num1) + mu0(num2) &
             & + gasr * temp * (61.781d0 - (10834.0d0 / temp)) &
             & - gasr * temp * dlog(1.0d12)
        
        property(s)%svap_cef%low%temp = temp
        property(s)%svap_cef%low%A = entr0_tmp
        property(s)%svap_cef%low%B = mu0_tmp + (entr0_tmp * temp) / 2.0d0 
        
        deallocate(entr0, mu0)        


        !上限を求める
        if (property(num1_cnd)%svap_cef%high%temp &
             & <= property(num2_cnd)%svap_cef%high%temp) then 
           temp = property(num1_cnd)%svap_cef%high%temp
        else
           temp = property(num2_cnd)%svap_cef%high%temp
        end if
        
        call entropy(temp, pres_ref, mol, property, entr0=entr0)
        entr0_tmp = entr0(num1) + entr0(num2) &
             & - gasr * 61.781d0  + gasr * dlog(1.0d12) 
        call entropy(temp - dt, pres_ref, mol, property, entr0=entr0_pre)
        entr0_tmp_pre = entr0_pre(num1) + entr0_pre(num2) &
             & - gasr * 61.781d0  + gasr * dlog(1.0d12) 
        call gibbs_energy(temp, pres_ref, mol, property, mu0=mu0)
        mu0_tmp  = mu0(num1) + mu0(num2) &
             & + gasr * temp * (61.781d0 - (10834.0d0 / temp)) &
             & - gasr * temp * dlog(1.0d12)
        
        entr0_dt = (entr0_tmp - entr0_tmp_pre)/ dt
        property(s)%svap_cef%high%temp = temp                
        property(s)%svap_cef%high%A = entr0_dt
        property(s)%svap_cef%high%B = entr0_tmp - entr0_dt * temp
        property(s)%svap_cef%high%C = mu0_tmp + entr0_tmp * temp &
             & - (entr0_dt * (temp ** 2.0d0)) / 2.0d0 
        
        deallocate(entr0, mu0, entr0_pre)
     end if
     
     deallocate(rvar)
     
  end do calc_svap

!!!
!!!後かたづけ
!!!
  deallocate(mol)

  
end subroutine def_svap_cef



