!---------------------------------------------------------------------
!     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_PressBasicZ(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)  :: ss_LatHeat(DimXMin:DimXMax, DimZMin:DimZMax) ! 
    real(8)  :: ss_MolWtV(DimXMin:DimXMax, DimZMin:DimZMax) !

    real(8)  :: Temp_0, Temp_1, Press_0, Press_1
                    ! 乾燥断熱線と湿潤断熱線が交わる高度を
                    ! 反復法を用いて計算する際に一時的に使う変数
    real(8)  :: LFC ! 乾燥断熱線と湿潤断熱線が交わる高度
    real(8)  :: LTP ! 湿潤断熱線と等温線が交わる高度
    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

       ss_TempAdia = TempSfc - Grav * ss_Z / ss_CpBasicZ
       ss_TempIso  = 135.0d0

       !--- 初期化. とりあえず乾燥断熱分布を与えておく
!       ss_TempBasicZ = ss_TempAdia
!       ss_ExnerBasicZ = ss_HydroStaticExnerZ(ss_TempBasicZ)
!!       ss_PressBasicZ = PressSfc * ss_ExnerBasicZ**(ss_CpBasicZ/GasR)
!       ss_PressBasicZ = 700.0d0
      
       !--- 下から圧力, 温度場を求め直す
!------(2006/01/13)-------------------------------------------
!
!       do k = RegZMin+1, DimZMax
!       do i = DimXMin, DimXMax
!
!         ss_TempSat(i,k) = SatPressB  &
!           & / (  &
!           &     SatPressA & 
!           &     - log(PressSfc & 
!           &           * ss_ExnerBasicZ(i,k)**(ss_CpBasicZ(i,k)/GasR) & 
!           &           / SatRatioCr &
!           &    )    &
!           & )  
!
!         if (ss_TempBasicZ(i,k) < ss_TempSat(i,k) ) then 
!  
!           do
!             ss_TempSat(i,k) = SatPressB  &
!               & / (  &
!               &     SatPressA & 
!               &     - log(  PressSfc & 
!               &           * ss_ExnerBasicZ(i,k)**(ss_CpBasicZ(i,k)/GasR) & 
!               &           / SatRatioCr) &
!               & )  
!
!             if ( abs(ss_TempSat(i,k) - ss_TempBasicZ(i,k)) > epsilon(1.0d0) ) then
!                ss_TempbasicZ(i,k) = ss_TempbasicZ(i,k) + (ss_TempSat(i,k) - ss_TempBasicZ(i,k))
!                ss_ExnerBasicZ(i,k) = ss_ExnerBasicZ(i,k-1) * &
!                     &  exp(- Grav / ss_CpBasicZ(i,k-1) * 0.5d0 * DelZ  &
!                     &  * (1.0d0 / ss_TempBasicZ(i,k-1) + 1.0d0 / ss_TempBasicZ(i,k)))
!             else
!               exit
!             end if
!           end do
!         end if
!
!         ss_TempBasicZ(i,k) = max(ss_TempBasicZ(i,k), ss_TempIso(i,k))
!      end do
!      end do
!------------------------------------------------

!------(2005/12/30)-----------------------------
!      do k = RegZMin+1, DimZMax
!         !--- 上の層の圧力を計算
!         ss_ExnerBasicZ(:,k) = ss_ExnerBasicZ(:,k-1)  &
!           & * exp(- Grav / ss_CpBasicZ(:,k) * DelZ  &
!           &       * 1.0d0 / ss_TempBasicZ(:,k-1) &
!           &      )
!         !--- 凝結温度を求める
!         !    ノイズによる凝結を防ぐため過飽和度を 1 % 小さくして計算
!         ss_TempSat(:,k) = SatPressB  &
!           & / (  &
!           &     SatPressA & 
!           &     - log(PressSfc & 
!           &           * ss_ExnerBasicZ(:,k)**(ss_CpBasicZ(:,k)/GasR) & 
!           &           / SatRatioCr &
!           &    )    &
!           &   )
!         !--- 上の層の温度を決める
!         ss_TempBasicZ(:,k) = &
!           & max(ss_TempAdia(:,k), ss_TempIso(:,k), ss_TempSat(:,k))
!       end do
!-----------------------------------------------

!------(2006/01/14)-----------------------------
!       ! (方法その 3) eccm を用いて温度分布を与える.
!
!       call eccm(1.0d0, ss_PressBasicZ(1,:), ss_TempBasicZ(1,:), ss_MolWtV(1,:), ss_LatHeat(1,:))
!
!       do k = DimZMin, DimZMax
!         ss_TempBasicZ(:,k) = ss_TempBasicZ(1,k)
!         ss_PressBasicZ(:,k) = ss_PressBasicZ(1,k)
!       end do
!       ss_ExnerBasicZ = (ss_PressBasicZ / PressSfc)**(GasR / ss_CpBasicZ)
!-----------------------------------

!-----------------------------------------------
!      ! (方法その 4) 乾燥断熱線, 湿潤断熱線, 等温線が交わる高度を計算し, 
!                     各領域で成り立つ式を用いて温度, 圧力を計算
!
       ! 乾燥断熱線と湿潤断熱線が交わる高度(LFC)を反復法で計算
       Press_0 = PressSfc
       Temp_0 = TempSfc
       do 
         Temp_1 = SatPressB / (SatPressA - dlog(Press_0/SatRatioCr))
         Press_1 = PressSfc*(Temp_1/TempSfc)**(ss_CpBasicZ(1,1)/GasR)
         if (abs(Temp_1 - Temp_0) < epsilon(0.0d0)) then
            LFC = TempSfc*ss_CpBasicZ(1,1)/Grav &
              &   * (1.0d0 - (Press_1/PressSfc)**(GasR/ss_CpBasicZ(1,1)))
           exit
         else
           Temp_0 = Temp_1
           Press_0 = Press_1
         end if
       end do

       ! 湿潤断熱線と唐音線が交わる高度(LTP)を計算
       LTP = LFC + GasR*SatPressB/Grav*dlog(Temp_1/ss_TempIso(1,1))

       ! LFC, LTP の値を表示. 
       ! write(*,*) 'LFC', LFC
       ! write(*,*) 'LTP', LTP

       do k = DimZMin, DimZMax
         if (s_Z(k) < LFC) then           ! 乾燥断熱線
           ss_TempBasicZ(:,k) = TempSfc - Grav/ss_CpBasicZ(:,k)*ss_Z(:,k)
           ss_PressBasicZ(:,k) =  &
             &  PressSfc*(ss_TempBasicZ(:,k)/TempSfc)**(ss_CpBasicZ(:,k)/GasR)
           ss_ExnerBasicZ(:,k) =  &
             &  ss_TempBasicZ(:,k)/TempSfc
         else if (s_Z(k) > LTP) then      ! 等温線
           ss_TempBasicZ(:,k) = ss_TempIso(:,k)
           ss_PressBasicZ(:,k) =   &
             &    exp(SatPressA - SatPressB/ss_TempIso(:,k))  &
             &  * exp(-Grav*(ss_Z(:,k) - LTP)/(GasR*ss_TempIso(:,k)))
           ss_ExnerBasicZ(:,k) =  &
             &    exp(GasR/ss_CpBasicZ(:,k)  &
             &        *(SatPressA - SatPressB/ss_TempIso(:,k)))  &
             &  * exp(-Grav*(ss_Z(:,k) - LTP)  &
             &        /(ss_CpBasicZ(:,k)*ss_TempIso(:,k))) &
             &  / PressSfc**(GasR/ss_CpBasicZ(:,k)) 
         else                             ! 湿潤断熱線
           ss_TempBasicZ(:,k) =  &
             &  Temp_1*exp(-Grav*(ss_Z(:,k) - LFC)/(GasR*SatPressB))
           ss_PressBasicZ(:,k) =   &
              &  SatRatioCr*exp(SatPressA - SatPressB/ss_TempBasicZ(:,k))
           ss_ExnerBasicZ(:,k) =   &
              &  (SatRatioCr/PressSfc)**(GasR/ss_CpBasicZ(:,k)) & 
              &  * exp(GasR/ss_CpBasicZ(:,k) &
              &        *(SatPressA - SatPressB/ss_TempBasicZ(:,k)))
         end if 
       end do

       ! 圧力から無次元圧力を計算
       ! ss_ExnerBasicZ = (ss_PressBasicZ/PressSfc)**(GasR/ss_CpBasicZ)

!--------------------------------------------------------

       !--- 温度, 圧力から温位を求める
       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

end module BasicSet