!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------

!= Module BasicSet
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: basicset.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!基本場の値を取得する. 
!   * basicset_file_init: 基本場の値を netCDF ファイルから取得
!   * basicset_calc_init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!   * 湿潤過程導入のため, 乾燥過程の時の基本場計算はもはや使わない. 
!
!== Future Plans
!

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

  !== save 属性
  save

  !== Public Interface
  real(8), allocatable  :: ss_DensBasicZ(:,:)     !密度
  real(8), allocatable  :: ss_PressBasicZ(:,:)    !無次元圧力
  real(8), allocatable  :: ss_ExnerBasicZ(:,:)    !無次元圧力
  real(8), allocatable  :: ss_TempBasicZ(:,:)     !温度
  real(8), allocatable  :: ss_PotTempBasicZ(:,:)  !温位
  real(8), allocatable  :: ss_VelSoundBasicZ(:,:) !音速
  real(8), allocatable  :: ss_MixRtVBasicZ(:,:)   !水蒸気混合比
  real(8), allocatable  :: ss_EffMolWtBasicZ(:,:) !分子量効果
  real(8), allocatable  :: ss_CpBasicZ(:,:)       !比熱
  

contains

  !== 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,  &! 配列の X 方向の下限
      &                   DimXMax,  &! 配列の X 方向の上限
      &                   DimZMin,  &! 配列の Z 方向の下限
      &                   DimZMax    ! 配列の Z 方向の上限

    
    !=== Input
    character(*), intent(in) :: basicfile
    
    character(20)  :: name            !変数名
    
    call BeginSub("basicset_file_init", &
      &                 fmt="%c",       &
      &                 c1="Get basic state variables from netCDF file.")
    
    !--- 配列の割り当て
    allocate(  &
      & ss_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),    &
      & ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),    &
      & ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),  &
      & ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_MixRtVBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),   & 
      & ss_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)        )
    
    !=== Get a Value from netCDF File
    name = "PressBasicZ"
    call HistoryGet( basicfile, name, ss_PressBasicZ )
    
    name = "ExnerBasicZ"
    call HistoryGet( basicfile, name, ss_ExnerBasicZ )

    name = "TempBasicZ"
    call HistoryGet( basicfile, name, ss_TempBasicZ )

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

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

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

    name = "MixRtVBasicZ"
    call HistoryGet( basicfile, name, ss_MixRtVBasicZ ) 

    name = "EffMolWtBasicZ"
    call HistoryGet( basicfile, name, ss_EffMolWtBasicZ ) 

    name = "CpBasicZ"
    call HistoryGet( basicfile, name, ss_CpBasicZ ) 

    call EndSub("basicset_file_init")

  end subroutine BasicSet_file_init


  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !典型的な基本場を計算するためのプログラム basic の初期値を
  !読み込み, 基本場の値を計算するためのモジュール
  !
  !湿潤断熱減率で温度, 分子量分布を与える. 
  !
  subroutine basicset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: DimXMin,  &! 配列の X 方向の下限
      &                   DimXMax,  &! 配列の X 方向の上限
      &                   DimZMin,  &! 配列の Z 方向の下限
      &                   DimZMax,  &! 配列の Z 方向の上限
      &                   RegZMin,  &! z 方向の物理領域の下限
      &                   s_Z,      &! Z 座標軸(スカラー格子点)
      &                   DelZ       ! 格子点間隔 [m]
    use physset,    only: PressSfc, &! 地表面圧力 [Pa]
      &                   TempSfc,  &! 地表面温度 [K]
      &                   Grav,     &! 重力加速度 [m/s^2]
      &                   MolWtDry, &! 乾燥成分の分子量 [kg/kg]
      &                   MolWtWet, &! 湿潤成分の分子量 [kg/kg]
      &                   GasRDry,  &! 乾燥成分の気体定数 [J/K kg]
      &                   CpDry,    &! 乾燥成分の平均定圧比熱 [J/K kg]
      &                   CvDry      ! 乾燥成分の平均定積比熱 [J/K kg]
    use bcset,    only: ss_BC
        
    !=== Input
    character(*), intent(in) :: cfgfile  !設定ファイル(NAMELIST)
    
    !=== Work
    integer  :: k
    real(8)  :: RH        ! 相対湿度
    real(8)  :: LCL       ! 凝結高度
    real(8)  :: TPL       ! 対流圏圏界面高度
    real(8)  :: Gamma     ! 温度減率
    real(8)  :: GammaWet  ! 温度減率
    real(8)  :: factor(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_SvapPressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)!飽和蒸気圧

    !=== NAMELIST 
    NAMELIST /basicset/ RH, LCL, GammaWet, TPL
    
    call BeginSub("basicset_calc_init", &
      &                 fmt="%c",       &
      &                 c1="calculate the basic state variables.")
    
    open (10, FILE=cfgfile)
    read(10, NML=basicset)
    close(10)
    
    !=== 変数の初期化
    allocate( &
      &    ss_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
      &    ss_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),    &
      &    ss_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),    &
      &    ss_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),     &
      &    ss_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),  &
      &    ss_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
      &    ss_MixRtVBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),   &
      &    ss_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), &
      &    ss_CpBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)        )

    !=== 初期化
    ss_DensBasicZ     = 0.0d0
    ss_PressBasicZ    = 0.0d0
    ss_ExnerBasicZ    = 0.0d0
    ss_TempBasicZ     = 0.0d0
    ss_PotTempBasicZ  = 0.0d0
    ss_VelSoundBasicZ = 0.0d0
    ss_MixRtVBasicZ   = 0.0d0
    ss_EffMolWtBasicZ = 0.0d0

    
    !=== 乾燥断熱温度分布, 静水圧平衡
    !==== 温度を計算. k=RegZMin が地表面. 対流圏界面より上は温度一定
    ss_TempBasicZ(:,RegZMin) = TempSfc
    Gamma = Grav / CpDry    

    do k = RegZMin, DimZMax-1      
      if (DelZ * k < LCL) then 
        ss_TempBasicZ(:,k+1) = ss_TempBasicZ(:,k) - Gamma * DelZ
      elseif (LCL <= DelZ * k .AND. DelZ * k < TPL) then 
        ss_TempBasicZ(:,k+1) = ss_TempBasicZ(:,k) - GammaWet * DelZ
      else
        ss_TempBasicZ(:,k+1) = ss_TempBasicZ(:,k) 
      end if
    end do
    !境界条件
    call boundary(ss_BC, ss_TempBasicZ)
!
!    do k = RegZMin, DimZMin+1, -1
!      ss_TempBasicZ(:,k-1) = ss_TempBasicZ(:,k) + Gamma * DelZ 
!    end do

    !==== 圧力を静水圧平衡から計算. k=0 が地表面
    ss_PressBasicZ(:,RegZMin) = PressSfc

    do k = RegZMin, DimZMax-1
      ss_PressBasicZ(:,k+1) = &
        & ss_PressBasicZ(:,k) &
        & - ss_PressBasicZ(:,k) * Grav * DelZ &
        &     / ( GasRDry * ss_TempBasicZ(:,k))
    end do
    !境界条件
    call boundary(ss_BC, ss_PressBasicZ)

!    do k = RegZMin, DimZMin+1, -1
!      ss_PressBasicZ(:,k-1) = &
!        & ss_PressBasicZ(:,k) &
!        & + ss_PressBasicZ(:,k) * Grav * DelZ &
!        &     / ( GasRDry * ss_TempBasicZ(:,k))
!    end do
    
    !=== 水蒸気混合比. 相対湿度 RH 一定. 
    !==== 飽和蒸気圧
    ss_SvapPressBasicZ = ss_svap_ss(ss_TempBasicZ)
        
    !=== 温位
    ss_PotTempBasicZ = &
      & ss_TempBasicZ * (PressSfc / ss_PressBasicZ) ** (GasRDry / CpDry) 

    !=== エクスナー関数
    ss_ExnerBasicZ = ss_TempBasicZ / ss_PotTempBasicZ    
    
    !==== 混合比
    !===== 相対湿度一定の場合    
    ss_MixRtVBasicZ = &
      & ((MolWtWet/MolWtDry) * ss_SvapPressBasicZ * RH) &
      & / (PressSfc * (ss_ExnerBasicZ **(CpDry / GasRDry)) &
      &                        - ss_SvapPressBasicZ * RH)
    
    !===== 相対湿度が単調減少の場合
!    write(*,*) "RH", RH
!
!    do k = DimZMin, DimZMax
!      factor(:,k) = RH - s_Z(k) * 1.0d-3 * 5.0d-2 
!      if (s_Z(k) < 0.0d0) factor(:,k)  = 0.0d0
!    end do
    
!    ss_MixRtVBasicZ = &
!      & ((MolWtWet/MolWtDry) * ss_SvapPressBasicZ * factor) &
!      & / (PressSfc * (ss_ExnerBasicZ ** (CpDry / GasRDry)) &
!      &                     - ss_SvapPressBasicZ * factor) 
!    where (ss_MixRtVBasicZ < 0.0d0) 
!      ss_MixRtVBasicZ = 0.0d0
!    end where
    
    !=== 分子量の効果
    ss_EffMolWtBasicZ = &
      &   (1.0d0 - (ss_MixRtVBasicZ / (MolWtWet/MolWtDry + ss_MixRtVBasicZ)))&
      & * (1.0d0 + ss_MixRtVBasicZ)
    
    !=== 密度
    ss_DensBasicZ = &
      & PressSfc * (ss_ExnerBasicZ ** (CvDry / GasRDry)) &
      &  / (GasRDry * ss_PotTempBasicZ / ss_EffMolWtBasicZ)
    
    !=== 音速
    ss_VelSoundBasicZ = &
      & sqrt ( &
      &   CpDry * GasRDry * ss_ExnerBasicZ * ss_PotTempBasicZ &
      &   / (CvDry * ss_EffMolWtBasicZ) &
      & )

    !=== 比熱. 他のサブルーチンとの整合性のため
    ss_CpBasicZ = CpDry
    
    call EndSub("basicset_calc_init")
    
  end subroutine basicset_calc_init


  function ss_svap_ss(ss_temp)
    use gridset,    only: DimXMin,  &! 配列の X 方向の下限
      &                   DimXMax,  &! 配列の X 方向の上限
      &                   DimZMin,  &! 配列の Z 方向の下限
      &                   DimZMax    ! 配列の Z 方向の上限

    ! 暗黙の型宣言禁止
    implicit none

    ! 入出力変数
    real(8)             :: ss_svap_ss(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !飽和蒸気圧 [Pa]
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !温度 [K]
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    ss_svap_ss = dexp( A - B / (C + ss_temp) ) !Antoine の式
    
  end function ss_svap_ss
  
end module BasicSet
