!= Module BasicEnvInit
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: initialdata_basic.f90,v 1.7 2011-06-03 06:11:21 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!デフォルトの基本場を設定するための変数参照型モジュール
!   * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得
!   * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
module initialdata_basic
  !
  ! デフォルトの基本場を作るためのルーチン

  !モジュール読み込み
  use dc_types,   only: STRING, DP
  use dc_message, only: MessageNotify
  use gridset,    only: kmin, kmax, &!配列の X 方向の上限・下限
    &                 ncmax             !凝縮成分の数
  use axesset, only:  z_Z,          &!スカラー格子点での高度
    &                 r_dz            !Z 方向の格子点間隔
  use constants, only: &
    &                 Grav,          &!重力加速度
    &                 TempSfc, PressSfc, &
    &                 GasRDry        !気体定数
  use composition, only: SpcWetMolFr     !凝縮成分の初期モル比
  use eccm,     only: ECCM_Dry,      &!
    &                 ECCM_Wet              
    
  !暗黙の型宣言禁止
  implicit none

  !デフォルトは private
  private

  real(8), save  :: Humidity = 0.0d0        !基本場の湿度 
!  real(8), save  :: TempSfc   = 300.0d0     !地表面温度       [K]
!  real(8), save  :: PressSfc  = 965.0d2     !地表面圧力       [Pa]
  real(8), save  :: TempStrat = 100.0d0     !成層圏の温度 [k]
  real(8), save  :: HeightStrat = 18.0d3    !成層圏の高度 [k]
  real(8), save  :: Dhight = 5.0d3          !重み関数のパラメータ [m]
  integer, save, public :: ID_FlaginitdataBasic = 0
  integer, parameter, public :: ID_FlaginitDataBasic_dry = 1
  integer, parameter, public :: ID_FlaginitDataBasic_moist = 2
  integer, parameter, public :: ID_FlaginitDataBasic_isothermal = 3

  ! public
  public initialdata_basic_nml
  public initialdata_basic_dry
  public initialdata_basic_moist
  public initialdata_basic_strat
  public initialdata_basic_isothermal
    
contains

  subroutine initialdata_basic_nml(cfgfile)
    ! nml の読み込み
    !
    use dc_iounit,     only : FileOpen    

    implicit none

    character(STRING), intent(in) :: cfgfile    
    integer                       :: unit     !装置番号
    character(STRING)  :: FlaginitdataBasic = ""
    
    NAMELIST /initialdata_basic/ &
      & FlaginitDataBasic, Humidity, TempStrat, HeightStrat, Dhight
    
    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=initialdata_basic)
    close(unit)

    if (FlaginitdataBasic == "moist") then 
      ID_FlaginitdataBasic = ID_FlaginitdataBasic_moist
    elseif (FlaginitdataBasic == "dry") then
      ID_FlaginitdataBasic = ID_FlaginitdataBasic_dry
    elseif (FlaginitdataBasic == "isothermal") then
      ID_FlaginitdataBasic = ID_FlaginitdataBasic_isothermal
    end if    

  end subroutine initialdata_basic_nml

  subroutine initialdata_basic_isothermal( z_TempBZ, z_PressBZ, za_MolFr )

    implicit none

    real(8), intent(out) :: z_TempBZ(kmin:kmax)
    real(8), intent(out) :: z_PressBZ(kmin:kmax)
    real(8), intent(out) :: za_MolFr(kmin:kmax, 1:ncmax)

    z_TempBZ = TempSfc
    z_PressBZ = PressSfc
    za_MolFr = 0.0d0

  end subroutine initialdata_basic_isothermal

  subroutine initialdata_basic_dry( z_TempBZ, z_PressBZ, za_MolFr )

    implicit none

    real(8), intent(out) :: z_TempBZ(kmin:kmax)
    real(8), intent(out) :: z_PressBZ(kmin:kmax)
    real(8), intent(out) :: za_MolFr(kmin:kmax, 1:ncmax)
    
    call eccm_dry( SpcWetMolFr(1:ncmax), Humidity, z_TempBZ, z_PressBZ, za_MolFr )

  end subroutine initialdata_basic_dry


  subroutine initialdata_basic_moist( z_TempBZ, z_PressBZ, za_MolFr )

    implicit none

    real(8), intent(out) :: z_TempBZ(kmin:kmax)
    real(8), intent(out) :: z_PressBZ(kmin:kmax)
    real(8), intent(out) :: za_MolFr(kmin:kmax, 1:ncmax)
    
    call eccm_wet( SpcWetMolFr(1:ncmax), Humidity, z_TempBZ, z_PressBZ, za_MolFr )
    
  end subroutine initialdata_basic_moist

  
  subroutine initialdata_basic_strat(z_TempBZ, z_PressBZ)
    !
    implicit none
    
    real(8), intent(inout) :: z_TempBZ(kmin:kmax)
    real(8), intent(inout) :: z_PressBZ(kmin:kmax)
    real(8) :: z_DTempDZ(kmin:kmax)
    real(8) :: DTempDZ
    real(8) :: Weight
    integer :: k    
    
    do k = 1, kmax-1
      z_DTempDZ(k) = (z_TempBZ(k) - z_TempBZ(k-1)) / r_dz(k-1)
    end do
    
    ! 対流圏界面より上の扱い
    !   圏界面より上は等温大気とする. 等温位大気から等温大気への遷移は 
    !   tanh を用いてなめらかにつなぐ
    do k = 1, kmax
      
      !重みつけの関数を用意. tanh を用いる
      Weight = ( tanh( (z_Z(k) - HeightStrat ) / Dhight ) + 1.0d0 ) * 5.0d-1
      
      !仮値として温度を計算する. 圏界面より上では TempStrat の等温大気に近づける
      z_TempBZ(k) = z_TempBZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
      
      !温度減率が断熱温度減率より小さくならないように
      DTempDZ = &
        & max( z_DTempDZ(k), (z_TempBZ(k) - z_TempBZ(k-1)) / r_dz(k-1) )
      
      !基本場の温度を決める
      z_TempBZ(k) = z_TempBZ(k-1) + DTempDZ * r_dz(k-1) 
      
      !圧力を静水圧平衡から計算
      z_PressBZ(k) =                                                    &
        &  z_PressBZ(k-1) * ( ( z_TempBZ(k-1) / z_TempBZ(k) )   &
        &    ** (Grav / ( DTempDZ * GasRDry ) ) )
    end do
    
  end subroutine initialdata_basic_strat
 
end module initialdata_basic
