!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module BasicSet
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: basicset.f90,v 1.2.2.6 2005/12/20 11:56:45 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

    !=== 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)  :: TempAdia ! 乾燥断熱線に沿った温度
    real(8)  :: TempSat  ! 凝結線に沿った温度
    real(8)  :: TempIso  ! 等温線に沿った温度
    integer  :: k
                                                                 !=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

! (2005/12/20 北守太一  書きかけ...)
!    elseif (TempType == "adia-sat-iso" .AND. PressType == "hydrostatic") then
!
!       !--- 最下層の圧力, 温度を計算
!       ss_ExnerBasicZ(:, RegZMin+1) = &
!         & exp(- Grav * s_Z(RegZMin+1) * 0.5d0   &
!         &     / (ss_CpBasicZ(:, RegZMin+1) * TempSfc))
!
!       ss_TempBasicZ(:, RegZMin+1) =  &
!         & TempSfc - Grav * s_Z(RegZMin+1) * 0.5d0 / ss_CpBasicZ(:, RegZMin+1)
!
!       TempAdia = TempSfc   &
!         &        - Grav * s_Z(RegZMin+1) * 0.5d0 / ss_CpBasicZ(1, RegZMin+1)
!       TempIso  = TempSfc - 25.0d0
!       TempSat  = SatPressB & 
!         & / (  &
!         &     (SatPressA - log(PressSfc)  &
!         &      - ss_CpBasicZ(1,RegZMin+1) / GasR  &
!         &        * log(ss_ExnerBasicZ(1,RegZMin+1)) &
!         &     )  &
!         &   ) 
!
!       !--- 下側から順番に圧力, 温度を計算
!      do k = RegZMin+2, DimZMax
!        ss_ExnerBasicZ(:,k) = ss_ExnerBasicZ(:,k)  &
!          & * exp(- Grav * DelZ / (ss_CpBasicZ(:,k-1) * ss_TempBasicZ(:,k-1)))
!
!         if (TempAdia > TempIso .and. TempAdia > TempSat) then
!           ss_TempBasicZ() = TempAdia
!
!         TempAdia = TempSfc   &
!           &        - Grav * s_Z(k) * 0.5d0 / ss_CpBasicZ(1,k)
!         TempIso  = TempSfc - 25.0d0
!         TempSat  = SatPressB & 
!           & / (  &
!           &     (SatPressA - log(PressSfc)  &
!           &      - ss_CpBasicZ(1,k) / GasR  &
!           &        * log(ss_ExnerBasicZ(1,k)) &
!           &     )  &
!           &   ) 
!       end do

    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.5d0 * DelZ  / TempSfc)

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

end module BasicSet

