!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module BasicSet
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: basicset.f90,v 1.2.2.8 2006/01/07 12:32:05 kitamo Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!基本場の値を取得する. 
!   * basicset_file_init: 基本場の値を netCDF ファイルから取得
!   * basicset_calc_init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
                                                                 !=end

module BasicSet
  !== 暗黙の型宣言禁止
  implicit none

  !== save 属性
  save
                                                                 !=begin
  !== Public Interface
  real(8), allocatable   :: ss_CpBasicZ(:,:)       !比熱
  real(8), allocatable   :: ss_DensBasicZ(:,:)     !密度
  real(8), allocatable   :: ss_ExnerBasicZ(:,:)    !無次元圧力
  real(8), allocatable   :: ss_PotTempBasicZ(:,:)    !温位
  real(8), allocatable   :: ss_VelSoundBasicZ(:,:) !音速
                                                                 !=end

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module
  !
  !基本場の値を netCDF ファイルから取得する. 
  !netCDF から値を取得するために gt4f90io を利用する. 
  !
  subroutine BasicSet_file_init(basicfile)

    !=== Dependency
    use dc_trace, only: BeginSub, EndSub
    use gt4_history
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax

    !=== Input
    character(*), intent(in) :: basicfile

                                                                 !=end
    
    character(20)  :: name            !変数名

    call BeginSub("basicset_file_init", &
&                 fmt="%c",       &
&                 c1="Get basic state variables from netCDF file.")
    
    !--- 配列の割り当て
    allocate(  &
         & ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) )

                                                                 !=begin
    !=== Get a Value from netCDF File
    name = "CpBasicZ"
    call HistoryGet( basicfile, name, ss_CpBasicZ )

    name = "ExnerBasicZ"
    call HistoryGet( basicfile, name, ss_ExnerBasicZ )

    name = "PotTempBasicZ"
    call HistoryGet( basicfile, name, ss_PotTempBasicZ )

    name = "DensBasicZ"
    call HistoryGet( basicfile, name, ss_DensBasicZ )

    name = "VelSoundBasicZ"
    call HistoryGet( basicfile, name, ss_VelSoundBasicZ ) 
                                                                 !=end   

    call EndSub("basicset_file_init")

  end subroutine BasicSet_file_init


                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !典型的な基本場を計算するためのプログラム basic の初期値を
  !読み込み, 基本場の値を計算するためのモジュール
  !
  !典型的な基本場のを計算するためのサブルーチン. 
  !タイプとしては以下が設定可能. 
  !   * 温度一様 or 温位一様
  !   * 圧力一様 or 静水圧平衡
  !
  subroutine basicset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: DimXMin, DimXMax, DimZMin, DimZMax, RegZMin, s_Z, DelZ
    use physset,    only: Grav, TempSfc, PressSfc, MolWtDry, CpDry, GasRUniv, GasR
    use cloudset,   only: SatPressA, SatPressB, SatRatioCr

    !=== Input
    character(*), intent(in) :: cfgfile  !設定ファイル(NAMELIST)
                                                                 !=end
    !=== Work
    character(80) :: TempType
    character(80) :: PressType
    real(8)  :: ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) ! 平均温度場
    real(8)  :: ss_CvBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_Z(DimXMin:DimXMax, DimZMin:DimZMax)          ! 2-D 座標
!!    real(8)  :: GasR
!    real(8)  :: s_TempAdiabat(DimZMin:DimZMax)            ! 等温時の温位分布
    real(8)  :: s_TempSaturate(DimZMin:DimZMax)            ! 等温時の温位分布
    real(8)  :: MidHeight                                       ! 等温と等温位が切り替わる高度(Odaka1998にて使用)
    real(8)  :: TempMidHeight                                   ! 高度 MidHeight での温度
    real(8)  :: ss_TempAdia(DimXMin:DimXMax, DimZMin:DimZMax) 
                                     ! 乾燥断熱線に沿った温度
    real(8)  :: ss_TempSat(DimXMin:DimXMax, DimZMin:DimZMax)
                                     ! 凝結線に沿った温度
    real(8)  :: ss_TempIso(DimXMin:DimXMax, DimZMin:DimZMax) 
                                     ! 等温線に沿った温度

    real(8)  :: s_TempBasicZ(DimZMin:DimZMax)
    real(8)  :: s_PresBasicZ(DimZMin:DimZMax)
    real(8)  :: s_MolWtV(DimZMin:DimZMax)
    real(8)  :: s_LatHeat(DimZMin:DimZMax)

    integer  :: k,i


                                                                 !=begin
    !=== NAMELIST 
    NAMELIST /basicset/ TempType, PressType
                                                                 !=end

    call BeginSub("basicset_calc_init", &
&                 fmt="%c",       &
&                 c1="calculate the basic state variables.")


    open (10, FILE=cfgfile)
    read(10, NML=basicset)
    close(10)

    !==== 確認
!    write(*,*) "TempType",  TempType
!    write(*,*) "PressType", PressType
    

    !=== 変数の初期化
    allocate( ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
         &    ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) )


    ss_CpBasicZ = 0.0d0
    ss_DensBasicZ = 0.0d0
    ss_ExnerBasicZ = 0.0d0
    ss_PotTempBasicZ = 0.0d0
    ss_VelSoundBasicZ = 0.0d0

    TempMidHeight = 0.0d0
    MidHeight = 0.0d0

    !--- 座標を初期化
    do k = DimZMin, DimZMax
       ss_Z(:,k) = s_Z(k)
    end do
    
    !--------------------------------------------
    !--- 気体定数を決める
    !--------------------------------------------
!!    GasR = GasRUniv / MolWtDry
    
    !--------------------------------------------
    !--- 比熱を決める
    !---   比熱は乾燥大気を念頭に一様とする. 
    !--------------------------------------------
    ss_CpBasicZ = CpDry / MolWtDry
    ss_CvBasicZ = ss_CpBasicZ - GasR
  
    !--------------------------------------------
    !--- 温度/温位, 無次元圧力の基本場
    !--------------------------------------------
    if (TempType == "uniform" .AND. PressType == "uniform") then 
       
       !--- 温度場は TempSfc で一様
       ss_TempBasicZ = TempSfc
       
       !--- 地表面圧力で一定. 定義より値は 1 となる. 
       ss_ExnerBasicZ = 1.0d0
       
       !--- 温位
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ 
       
    elseif (TempType == "uniform" .AND. PressType == "hydrostatic") then 
       
       !--- 温度場は TempSfc で一様
       ss_TempBasicZ = TempSfc
       
       !--- 無次元圧力と温位を連立して解く. 
       ss_ExnerBasicZ = &
            & dexp( - Grav * ss_Z / ( ss_CpBasicZ * ss_TempBasicZ ) )

       !--- 無次元圧力を利用して温位を計算
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ
       
    elseif (TempType == "adiabat" .AND. PressType == "uniform") then 
       
       !--- 地表面での温位(= TempSfc)で一定
       ss_PotTempBasicZ = TempSfc
       
       !--- 地表面圧力で一定. 定義より値は 1 となる. 
       ss_ExnerBasicZ = 1.0d0
       
    elseif (TempType == "adiabat" .AND. PressType == "hydrostatic") then 
       
       !--- 地表面での温位(= TempSfc)で一定
       ss_PotTempBasicZ = TempSfc
       
       !--- 無次元圧力を計算. 温位一定
       ss_ExnerBasicZ = &
            & - Grav * ss_Z / ( ss_CpBasicZ * ss_PotTempBasicZ ) + 1.0d0

    elseif (TempType == "odaka1998" .AND. PressType == "hydrostatic") then 

       !--- 温度分布を与える. 各高度毎に等温線と乾燥断熱線の高い方を選ぶ.
       ss_TempBasicZ =   &
         & max(TempSfc - 25.0d0, TempSfc - Grav * ss_Z / ss_CpBasicZ) 

       !--- 与えた温度分布から静水圧平衡の式を用いて圧力を求める
       ss_ExnerBasicZ = ss_HydroStaticExnerZ(ss_TempBasicZ)

       !--- 温度, 圧力から温位を求める
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ

    elseif (TempType == "adia-sat-iso" .AND. PressType == "hydrostatic") then

       !--- 温度, 圧力分布を求める
       call eccm(1.0d0, s_PresBasicZ, s_TempBasicZ, s_MolWtV, s_LatHeat)

       do k = DimZMin, DimZMax
         ss_TempBasicZ(:,k) = s_TempBasicZ(k)
         ss_ExnerBasicZ(:,k) =  (s_PresBasicZ(k) / PressSfc)**(GasR/ss_CpBasicZ(RegZMin,k))
       end do
       !--- 温度, 圧力から温位を求める
       ss_PotTempBasicZ = ss_TempBasicZ / ss_ExnerBasicZ

    else
       call MessageNotify("Message", "basicset_calc_init", &
            &             "unknown TempType, or PressType")
    end if
    
    
    !--------------------------------------------
    !---- 密度
    !--------------------------------------------
    ss_DensBasicZ = &
         & PressSfc * ( ss_ExnerBasicZ ** (ss_CvBasicZ / GasR) ) &
         &   / ( GasR * ss_PotTempBasicZ )
    
    
    !--------------------------------------------    
    !--- 平均場の音速
    !--------------------------------------------
    ss_VelSoundBasicZ = &
         & ss_CpBasicZ * GasR * ss_ExnerBasicZ * ss_PotTempBasicZ &
         &   / ss_CvBasicZ  
    ss_VelSoundBasicZ = sqrt( ss_VelSoundBasicZ )


    !==== 確認
    ! write(*,*) "ss_CpBasicZ ", maxval( ss_CpBasicZ )
    ! write(*,*) "ss_DensBasicZ ", maxval( ss_DensBasicZ )
    ! write(*,*) "ss_ExnerBasicZ ", maxval( ss_ExnerBasicZ )
    ! write(*,*) "ss_PotTempBasicZ ", maxval( ss_PotTempBasicZ )
    ! write(*,*) "ss_VelSoundBasicZ ", maxval( ss_VelSoundBasicZ )

    call EndSub("basicset_calc_init")

  end subroutine basicset_calc_init

  function ss_HydroStaticExnerZ(ss_TempZ)
                                                                  !=begin
  !== Overviews
  ! 静水圧平衡を仮定して, 温度分布から無次元圧力分布を求める. 
  ! 地表面から上端まで静水圧平衡の式を鉛直積分することで計算する. 
  !
  !== Dependancy 
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, RegZMin, DelZ
    use physset, only: PressSfc, TempSfc, Grav

  !==Input
    real(8), intent(in)     :: ss_TempZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 ! 温度

  !==Output
    real(8)                 :: ss_HydroStaticExnerZ(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 ! 無次元圧力
                                                                  !=end
  !==Work
    integer                 :: k

    ! 初期化. 地表面より下の圧力は計算しないのでこのまま 1 が入る. 
    ss_HydroStaticExnerZ = 1.0d0    

    ! 最下層の圧力を計算
    ss_HydroStaticExnerZ(:,RegZMin+1) = &
      &   exp(- Grav / ss_CpBasicZ(:,RegZMin+1) * 0.25d0 * DelZ     &
      &       * (1.0d0 / TempSfc + 1.0d0 / ss_TempZ(:,RegZMin+1) )  &
      &      )

    ! ひとつ下の格子点の圧力から次の格子点の圧力を求める
    do k = RegZMin+2, DimZMax
      ss_HydroStaticExnerZ(:,k) = ss_HydroStaticExnerZ(:,k-1) &
        & * exp(- Grav / ss_CpBasicZ(:,k) * 0.5d0 * DelZ  &
        &       * (1.0d0 / ss_TempZ(:,k-1) + 1.0d0 / ss_TempZ(:,k)) &
        &      )
    end do

    ! 下側マージン部分の圧力を計算
    do k = RegZMin, DimZMin, -1
      ss_HydroStaticExnerZ(:,k) = ss_HydroStaticExnerZ(:,k+1) &
        & * exp(Grav / ss_CpBasicZ(:,k) * 0.5d0 * DelZ  &
        &       * (1.0d0 / ss_TempZ(:,k+1) + 1.0d0 / ss_TempZ(:,k)) &
        &      )
    end do
  end function ss_HydroStaticExnerZ

!= Subroutine ECCM
!
!   * Developer: SUGIYAMA Ko-ichiro 
!   * Version: $Id: eccm.f90,v 1.1.2.1 2006/01/13 06:10:20 sugiyama Exp $
!   * Tag Name: $Name: arare3j $
!   * Change History: 
!
!== Overview 
!
!断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!

subroutine eccm( MolFrIni, Pres, Temp, MolFrV, LatentHeat )
  use gridset, only: DimZMin, DimZMax, RegZMin, DelZ
  use physset, only: GasRUniv, TempSfc, PressSfc, Grav, &
                     & MolWtDry, CpDry ! CpWetMol, MolWtWet

  implicit none

  real(8), intent(in)      :: MolFrIni
  real(8), intent(out)     :: Temp( DimZMin:DimZMax )
  real(8), intent(out)     :: Pres( DimZMin:DimZMax )
!  real(8), intent(out)     :: MixRat( DimZMin:DimZMax )
  real(8), intent(out)     :: MolFrV( DimZMin:DimZMax )
  real(8), intent(out)     :: LatentHeat( DimZMin:DimZMax )

!!!
!!!設定する必要のある変数
!!!
  ! H2O(g,l) の値
  real(8), parameter       :: Antoine_A = 27.3d0
  real(8), parameter       :: Antoine_B = 3103.0d0
  real(8), parameter       :: Antoine_C = 0.0d0

!!!
!!!設定する必要がない変数
!!!
  real(8)                  :: MolFr( DimZMin:DimZMax, 1:2 ) !モル分率
  real(8)                  :: Cp( 1:2 )                     !比熱   (J/K mol)
  real(8)                  :: CpMean( DimZMin:DimZMax )     !比熱   (J/K mol)
  real(8)                  :: MolWtMean( DimZMin:DimZMax )  !分子量 
  real(8)                  :: MolWt( 1:2 )                  !分子量 
  
  real(8)                  :: VapPres           !分圧
  real(8)                  :: SatVapPres        !飽和蒸気圧
  
  real(8)                  :: dtdz
  real(8)                  :: svap
  real(8)                  :: A, B
  
  integer                  :: i
  integer                  :: d = 1
  integer                  :: w = 2
  real(8)                  :: MolWtWet
  real(8)                  :: CpDryMol
  real(8)                  :: CpWetMol

  MolWtWet = MolWtDry
  CpDryMol = CpDry
  CpWetMol = CpDry

!!!
!!!配列の初期化
!!!
  Temp = TempSfc
  MolFr = 0.0d0
  CpMean = 0.0d0 
  LatentHeat = 0.0d0

  !比熱
  Cp(d) = CpDryMol
  Cp(w) = CpWetMol

  !分子量
  MolWt(d) = MolWtDry
  MolWt(w) = MolWtWet

  !初期モル比
  MolFr(:,w) = MolFrIni 
  MolFr(:,d) = 1.0d0 - MolFr(:,w) 

    !地表面での温度・圧力
    Pres = 0.0d0
    Temp = 0.0d0
    Pres( RegZMin ) = PressSfc
    Temp( RegZMin ) = TempSfc

  !!!
  !!!断熱減率 dT/dz の計算. 
  !!!
    do i = RegZMin, DimZMax-1

       !前のステップでのモル分率を用いて現在の蒸気圧を計算
       ! VapPres = MolFr(i-1, w) * Pres(i)
       VapPres = Pres(i)

       !飽和蒸気圧
       call antoine( Temp(i), svap = SatVapPres )

       !飽和蒸気圧から凝結の有無を決める
!       if ( VapPres < SatVapPres ) then 

          !モル比
          MolFr(i,:) = MolFr(i-1,:)

          !平均分子量の計算 (モル数は前のステップと変わらない)
          MolWtMean(i) = dot_product( MolFr(i, d:w), MolWt(d:w) )

          !平均比熱の計算 (モル数は前のステップと変わらない)
          CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )

          !乾燥断熱減率を計算
          dtdz = - MolWtMean(i) * Grav / CpMean(i)

!       else
!
!          !飽和蒸気圧
!          call Antoine( Temp(i), svap )
!
!          !潜熱
!          call PhaseChangeEntalpy(Temp(i), LatentHeat(i))
!
!          !飽和蒸気圧と圧力から現在のモル数を計算
!          MolFr(i,w) = svap / Pres(i)
!          MolFr(i,d) = 1.0d0 - MolFr(i,w)
!
!          !平均比熱の計算
!          CpMean(i)    = dot_product( MolFr(i,d:w), Cp(d:w) )
!
!          !平均分子量の計算
!          MolWtMean(i) = dot_product( MolFr(i,d:w), MolWt(d:w) )
!
!          !係数. 温度 Temp(i) で評価
!          A = LatentHeat(i) * MolFr(i,w) / ( GasRUniv * Temp(i) )
!          B = ( LatentHeat(i) ** 2.0d0 ) * MolFr(i,w)  &
!               & / ( CpMean(i) * GasRUniv * ( Temp(i) ** 2.0d0 ) )
!
!          !断熱温度減率
!          dtdz = - MolWtMean(i) * Grav * ( 1.0d0 + A ) / ( CpMean(i) * ( 1.0d0 + B ) )
!
!       end if

       !温度を計算
       Temp(i+1) = Temp(i) + dtdz * DelZ
       write(*,*) Temp(i+1)

       !圧力を静水圧平衡から計算
       Pres(i+1) = &
            & Pres(i)  - ( Pres(i) * Grav * DelZ * MolWtDry ) / ( GasRUniv * 0.5d0* (Temp(i)+Temp(i+1)) )

       !モル比 (DimZMax でも値を持つようにするため)
       MolFr(i+1,:) = MolFr(i,:)

    end do

    do i = RegZMin - 1, DimZMin, - 1 

       !モル比
       MolFr(i,:) = MolFr(i+1,:)

       !平均分子量の計算 (モル数は前のステップと変わらない)
       MolWtMean(i) = dot_product( MolFr(i+1, d:w), MolWt(d:w) )

       !平均比熱の計算 (モル数は前のステップと変わらない)
       CpMean(i) = dot_product( MolFr(i+1, d:w), Cp(d:w) )

       !乾燥断熱減率を計算
       dtdz = MolWtMean(i+1) * Grav / CpMean(i)

       !温度を計算
       Temp(i) = Temp(i+1) + dtdz * DelZ

       !圧力を静水圧平衡から計算
       Pres(i) = &
            & Pres(i+1)  + ( Pres(i+1) * Grav * DelZ * MolWtDry ) / ( GasRUniv * Temp(i+1) )

    end do

    !=== 出力用
    MolFrV = MolFr(:,w)
    LatentHeat(:) = LatentHeat(:) / MolWt(w)

  end subroutine eccm

  subroutine Antoine(temp, svap)    
    implicit none
    real(8), parameter       :: Antoine_A = 27.3d0
    real(8), parameter       :: Antoine_B = 3103.0d0
    real(8), parameter       :: Antoine_C = 0.0d0
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: svap         !飽和蒸気圧
    real(8)                :: svap_log
  
    !Antoine の式
!    svap_log = (Antoine_A - (Antoine_B / (Antoine_C + Temp - 273.15D0) ) ) &
!         & * dlog(10.0D0) + dlog(133.322D0)
    svap_log = (Antoine_A - (Antoine_B / (Antoine_C + Temp) ) )
    svap = dexp(svap_log)
    
  end subroutine Antoine


  subroutine PhaseChangeEntalpy(temp, LatentHeat)    
    use physset, only: GasRUniv
    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: LatentHeat   !潜熱
    real(8)                :: svap_log_dt
    real(8), parameter       :: Antoine_A = 27.3d0  
    real(8), parameter       :: Antoine_B = 3103.0d0
    real(8), parameter       :: Antoine_C = 0.0d0   

!    svap_log_dt =  antoine_B * dlog(10.0D0) &
!         & / (antoine_C + temp - 273.15D)**2.0D0 
    svap_log_dt =  antoine_B * dlog(10.0D0) &
         & / (antoine_C + temp)**2.0D0 
    
    LatentHeat = svap_log_dt * GasRUniv * ( Temp ** 2.0d0 ) 
 
  end subroutine PhaseChangeEntalpy
  
end module BasicSet