!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module BasicSet
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: basicset.f90,v 1.2.2.3 2005/09/16 14:48:13 kitamo Exp $
!   * Tag Name: $Name: arare3m $
!   * 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, s_Z
    use physset,    only: Grav, TempSfc, PressSfc, MolWtDry, CpDry, GasRUniv, GasR
    
    !=== 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_PotTempIsothermal(DimZMin:DimZMax)            ! 等温時の温位分布
    real(8)  :: MidHeight                                       ! 等温と等温位が切り替わる高度(Odaka1998にて使用)
    real(8)  :: TempMidHeight                                   ! 高度 MidHeight での温度
    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 

       !---高度 5 km での温度を計算
       MidHeight = 5000.0d0   
       TempMidHeight = TempSfc - Grav * MidHeight / ss_CpBasicZ(DimXMax, DimZMax) 

       s_PotTempIsothermal(:) =     & 
             & TempMidHeight * exp( Grav * s_Z(:) / ( ss_CpBasicZ(DimXMax, :) * TempMidHeight ) )

       do k = DimZMin, DimZMax
          if ( TempSfc > s_PotTempIsothermal(k) ) then
             ss_PotTempBasicZ(:,k) = TempSfc
             ss_ExnerBasicZ(:,k) = &
                  & - Grav * ss_Z(:,k) / ( ss_CpBasicZ(:,k) * ss_PotTempBasicZ(:,k) ) + 1.0d0
          else
             ss_PotTempBasicZ(:,k) = s_PotTempIsothermal(k)
             ss_ExnerBasicZ(:,k) = &
                  & - Grav * ss_Z(:,k) / ( ss_CpBasicZ(:,k) * ss_PotTempBasicZ(:,k) ) + 1.0d0
          end if
             
!          write(*,*) s_Z(k), ss_PotTempBasicZ(DimXMax,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
  
end module BasicSet
