!= Module BasicEnvInit
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: initialdata_basic.f90,v 1.1 2011-04-26 04:08:10 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_3d,  only: kmin, kmax, &!配列の X 方向の上限・下限
    &                 na,            &!凝縮成分の数
    &                 z_Z,           &!スカラー格子点での高度
    &                 r_dz            !Z 方向の格子点間隔
  use basicset_3d, only: &
    &                 MolWtDry,      &!乾燥成分の分子量
    &                 Grav,          &!重力加速度
    &                 SpcWetMolFr,   &!凝縮成分の初期モル比
    &                 GasRUniv        !気体定数
  use eccm_3d,  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_FlagInitialdataBasic = 0
  integer, parameter, public :: ID_FlagInitialDataBasic_dry = 1
  integer, parameter, public :: ID_FlagInitialDataBasic_moist = 2

  ! public
  public initialdata_basic_nml
  public initialdata_basic_dry
  public initialdata_basic_moist
  public initialdata_basic_strat
    
contains

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

    implicit none

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

    if (FlagInitialdataBasic == "moist") then 
      ID_FlagInitialdataBasic = ID_FlagInitialdataBasic_moist
    elseif (FlagInitialdataBasic == "dry") then
      ID_FlagInitialdataBasic = ID_FlagInitialdataBasic_dry
    end if    

  end subroutine initialdata_basic_nml


  subroutine initialdata_basic_dry( z_TempBasicZ, z_PressBasicZ, za_MolFr )

    implicit none

    real(8), intent(out) :: z_TempBasicZ(kmin:kmax)
    real(8), intent(out) :: z_PressBasicZ(kmin:kmax)
    real(8), intent(out) :: za_MolFr(kmin:kmax, 1:na)
    
    call eccm_dry( SpcWetMolFr(1:na), Humidity, z_TempBasicZ, z_PressBasicZ, za_MolFr )

  end subroutine initialdata_basic_dry


  subroutine initialdata_basic_moist( z_TempBasicZ, z_PressBasicZ, za_MolFr )

    implicit none

    real(8), intent(out) :: z_TempBasicZ(kmin:kmax)
    real(8), intent(out) :: z_PressBasicZ(kmin:kmax)
    real(8), intent(out) :: za_MolFr(kmin:kmax, 1:na)
    
    call eccm_wet( SpcWetMolFr(1:na), Humidity, z_TempBasicZ, z_PressBasicZ, za_MolFr )
    
  end subroutine initialdata_basic_moist

  
  subroutine initialdata_basic_strat(z_TempBasicZ, z_PressBasicZ)
    !
    implicit none
    
    real(8), intent(inout) :: z_TempBasicZ(kmin:kmax)
    real(8), intent(inout) :: z_PressBasicZ(kmin:kmax)
    real(8) :: z_DTempDZ(kmin:kmax)
    real(8) :: DTempDZ
    real(8) :: Weight
    real(8) :: tmp
    integer :: k    
    
    do k = 1, kmax-1
      z_DTempDZ(k) = (z_TempBasicZ(k) - z_TempBasicZ(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 の等温大気に近づける
      tmp =  z_TempBasicZ(k) 
      z_TempBasicZ(k) = z_TempBasicZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
!      write(*,*) Weight, tmp, "=>",  z_TempBasicZ(k)
      
      !温度減率が断熱温度減率より小さくならないように
      DTempDZ = &
        & max( z_DTempDZ(k), (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / r_dz(k-1) )
      
      !基本場の温度を決める
      z_TempBasicZ(k) = z_TempBasicZ(k-1) + DTempDZ * r_dz(k-1) 
      
      !圧力を静水圧平衡から計算
      z_PressBasicZ(k) =                                                    &
        &  z_PressBasicZ(k-1) * ( ( z_TempBasicZ(k-1) / z_TempBasicZ(k) )   &
        &    ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
    end do
    
  end subroutine initialdata_basic_strat
  

end module initialdata_basic

