module saturate

contains

  subroutine xyz_CalcQVapSat( Temp1, Temp2, MixPhaseType, xyz_Temp, xyz_Press, FillValue, xyz_QVapSat )
    !
    ! 温度 *Temp* と気圧 *Press* を用い, 
    ! 飽和比湿 *QVapSat* を求めます. 
    !
    ! Calculate saturation specific humidity *QVapSat* using
    ! temperature *Temp* and air pressure *Press*. 
    !

    ! モジュール引用 ; USE statements
    !
    use vtype_module


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: Temp1
    real(DP), intent(in):: Temp2
    character(*), intent(in) :: MixPhaseType
    real(DP), intent(in):: xyz_Temp (:,:,:)
                              ! $ T $ . 温度. Temperature
    real(DP), intent(in):: xyz_Press(:,:,:)
                              ! $ p $ . 気圧. Air pressure
    real(DP), intent(in):: FillValue

    real(DP), intent(out):: xyz_QVapSat(:,:,:)
                              ! $ q^{*} $ . 飽和比湿. Saturation specific humidity


    ! 作業変数
    ! Work variables
    !
    real(dp), parameter:: Es0 = 611.0_DP
                              ! 0 ℃の飽和蒸気圧 [Pa]. 
                              ! Saturation water vapor pressure at 0 deg C [Pa]
    real(dp), parameter:: GasRUniv = 8.314_DP
    real(dp), parameter:: MolWtDry = 28.964e-3_DP
    real(dp), parameter:: MolWtWet = 18.01528e-3_DP
    real(dp), parameter:: LatentHeat = 2.5e6_DP
    real(dp), parameter:: LatentHeatFusion = 334.0e3_DP
!!$    real(dp), parameter:: LatentHeat = 2.5e6_DP + 334.0d3

    real(dp) :: GasRWet
    real(dp) :: EpsV

    real(DP) :: xyz_WatFrac(size(xyz_Temp,1),size(xyz_Temp,2),size(xyz_Temp,3))
    real(DP) :: xyz_QVapSatOnLiq(size(xyz_Temp,1),size(xyz_Temp,2),size(xyz_Temp,3))
    real(DP) :: xyz_QVapSatOnSol(size(xyz_Temp,1),size(xyz_Temp,2),size(xyz_Temp,3))

    integer :: i, j, k

    GasRWet          = GasRUniv / MolWtWet
    EpsV             = MolWtWet / MolWtDry

    ! 実行文 ; Executable statement
    !

    call CalcWatFrac( Temp1, Temp2, MixPhaseType, xyz_Temp, xyz_WatFrac )

    do k = 1, size( xyz_Temp, 3 )
      do j = 1, size( xyz_Temp, 2 )
        do i = 1, size( xyz_Temp, 1 )
          if ( xyz_Temp(i,j,k) /= FillValue ) then
            xyz_QVapSatOnLiq(i,j,k) = &
              & EpsV * Es0 * exp( LatentHeat / GasRWet * ( 1./273. - 1./xyz_Temp(i,j,k) ) ) &
              & / xyz_Press(i,j,k)
            xyz_QVapSatOnSol(i,j,k) = &
              & EpsV * Es0 * exp( ( LatentHeat + LatentHeatFusion ) / GasRWet * ( 1./273. - 1./xyz_Temp(i,j,k) ) ) &
              & / xyz_Press(i,j,k)
            xyz_QVapSat(i,j,k) = xyz_WatFrac(i,j,k) * xyz_QVapSatOnLiq(i,j,k) + ( 1.0_DP - xyz_WatFrac(i,j,k) ) * xyz_QVapSatOnSol(i,j,k)
          else
            xyz_QVapSat(i,j,k) = FillValue
          end if
        end do
      end do
    end do


  end subroutine xyz_CalcQVapSat

  subroutine CalcWatFrac( Temp1, Temp2, MixPhaseType, xyz_Temp, xyz_WatFrac )

    ! モジュール引用 ; USE statements
    !
    use vtype_module


    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: Temp1
    real(DP), intent(in):: Temp2
    character(*), intent(in) :: MixPhaseType
    real(DP), intent(in):: xyz_Temp (:,:,:)
    real(DP), intent(out):: xyz_WatFrac(:,:,:)


    ! 作業変数
    ! Work variables
    !

    integer :: i, j, k

    ! 実行文 ; Executable statement
    !

    if ( MixPhaseType == 'Lin' ) then

      if ( Temp1 == Temp2 ) then
        do k = 1, size( xyz_Temp, 3 )
          do j = 1, size( xyz_Temp, 2 )
            do i = 1, size( xyz_Temp, 1 )
              if ( xyz_Temp(i,j,k) > Temp1 ) then
                xyz_WatFrac(i,j,k) = 1.0_DP
              else
                xyz_WatFrac(i,j,k) = 0.0_DP
              end if
            end do
          end do
        end do
      else
        xyz_WatFrac = ( 1.0_DP - 0.0_DP ) / ( Temp2 - Temp1 ) * ( xyz_Temp - Temp1 ) + 0.0_DP
        xyz_WatFrac = min( max( xyz_WatFrac, 0.0_DP ), 1.0_DP )
      end if
    else if ( MixPhaseType == 'Quad' ) then
      xyz_WatFrac = ( max( ( xyz_Temp - Temp1 ), 0.0_DP ) / ( Temp2 - Temp1 ) )**2
      xyz_WatFrac = min( max( xyz_WatFrac, 0.0_DP ), 1.0_DP )
    else
      write( 6, * ) 'Unexpected type: ', trim( MixPhaseType )
      stop
    end if


  end subroutine CalcWatFrac

end module saturate
