!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine BasicEnv
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!デフォルトの基本場を設定するための変数参照型モジュール
!   * BasicEnvFile_init: 基本場の値を netCDF ファイルから取得
!   * BasicEnvCalc_Init: 基本場の情報を Namelist から取得して値を計算
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!

subroutine BasicEnv()
  !
  !デフォルトの基本場を設定するためのサブルーチン. 
  !基本場を計算し, BasicSet モジュールの値を初期化する. 
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している. 
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである. 
  !

  !モジュール読み込み
  use gridset,  only: DimXMin,       &!配列の X 方向の下限
    &                 DimXMax,       &!配列の X 方向の上限
    &                 DimZMin,       &!配列の Z 方向の下限
    &                 DimZMax,       &!配列の Z 方向の上限
    &                 SpcNum,        &!凝縮成分の数
    &                 s_Z,           &!スカラー格子点での高度
    &                 DelZ            !Z 方向の格子点間隔
  use basicset, only: BasicSetArray_Init, &!
    &                 PressSfc,      &!乾燥成分の定圧比熱  
    &                 GasRDry,       &!乾燥成分の定圧比熱
    &                 CpDry,         &!乾燥成分の定圧比熱
    &                 CvDry,         &!乾燥成分の定積比熱
    &                 MolWtDry,      &!乾燥成分の分子量
    &                 Grav,          &!重力加速度
    &                 SpcWetMolFr,   &!凝縮成分の初期モル比
    &                 MolWtWet,      &!凝縮成分の分子量
    &                 EnvType,       &!基本場の設定
    &                 Tropopause,    &!対流圏界面高度
    &                 Humidity        !基本場の相対湿度
  use Boundary, only: xz_Boundary_Cyc_Sym 
  use MoistFunc,only: xz_MolFr2MixRt
  use ECCM,     only: ECCM_Init,     &!
    &                 ECCM_Temp_Press, &!
    &                 ECCM_MolFr
  
  !暗黙の型宣言禁止
  implicit none

  !変数の定義
  real(8) :: xz_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_PressBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_TempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_VelSoundBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  real(8) :: xz_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)

  real(8) :: MolFrIni(SpcNum)
  real(8) :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  real(8) :: xz_MixRtDivMolWt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  integer :: i, k, s

  !---------------------------------------------------------------
  ! 配列の初期化
  !---------------------------------------------------------------
  xz_DensBasicZ     = 0.0d0
  xz_PressBasicZ    = 0.0d0
  xz_ExnerBasicZ    = 0.0d0
  xz_TempBasicZ     = 0.0d0
  xz_PotTempBasicZ  = 0.0d0
  xz_VelSoundBasicZ = 0.0d0
  xz_MixRtBasicZ    = 0.0d0
  xz_EffMolWtBasicZ = 0.0d0

  !---------------------------------------------------------------
  ! ECCM モジュールの初期化
  !---------------------------------------------------------------
  call ECCM_Init
       
  !---------------------------------------------------------------
  ! EnvType を元に, 温度, 圧力, 組成を決める
  !---------------------------------------------------------------
  if (trim(EnvType) == 'Dry') then 
    MolFrIni = 0.0d0
  elseif (trim(EnvType) == 'Moist') then 
    MolFrIni = SpcWetMolFr(1:SpcNum) 
  end if
 
  do i = DimXMin, DimXMax
    call ECCM_Temp_Press( MolFrIni, xz_TempBasicZ(i,:), xz_PressBasicZ(i,:) ) 
  end do
  
  ! 対流圏界面より上の扱い
  do k = DimZMin, DimZMax      
    if (s_Z(k) > Tropopause) then 
      xz_TempBasicZ(:,k) = xz_TempBasicZ(:,k-1) 
      xz_PressBasicZ(:,k) = &
        & xz_PressBasicZ(:,k-1) &
        &   * dexp( - DelZ * Grav / (GasRDry * xz_TempBasicZ(:,k)) ) 
    end if
  end do
  
  !境界条件
  xz_TempBasicZ = xz_Boundary_Cyc_Sym(xz_TempBasicZ)
  xz_PressBasicZ = xz_Boundary_Cyc_Sym(xz_PressBasicZ)
  
  !---------------------------------------------------------------
  ! 混合比
  !---------------------------------------------------------------
  do i = DimXMin, DimXMax      
    call ECCM_MolFr( SpcWetMolFr(1:SpcNum) * Humidity, Humidity, xz_TempBasicZ(i,:), &
      &              xz_PressBasicZ(i,:), xz_MolFr(i,:,:) )
  end do
    
  !気相のモル比を混合比に変換
  do s = 1, SpcNum
    xz_MixRtBasicZ(:,:,s) = xz_MolFr2MixRt(MolWtWet(s), xz_MolFr(:,:,s))
  end do

  !値が小さくなりすぎないように最低値を与える
  where (xz_MixRtBasicZ <= 1.0d-20 )
    xz_MixRtBasicZ = 1.0d-20
  end where
  
  !境界条件
  do s = 1, SpcNum
    xz_MixRtBasicZ(:,:,s) = xz_Boundary_Cyc_Sym(xz_MixRtBasicZ(:,:,s))
  end do

      
  !---------------------------------------------------------------
  ! 分子量の効果
  !---------------------------------------------------------------
  do s = 1, SpcNum
    xz_MixRtDivMolWt(:,:,s) = xz_MixRtBasicZ(:,:,s) / MolWtWet(s)
  end do
  
  xz_EffMolWtBasicZ = &
    & 1.0d0 -  sum(xz_MixRtDivMolWt,3) &
    &            / ((1.0d0 / MolWtDry) + sum(xz_MixRtDivMolWt,3))  &
    & * (1.0d0 + sum(xz_MixRtBasicZ,3) )
  
  xz_EffMolWtBasicZ = xz_Boundary_Cyc_Sym(xz_EffMolWtBasicZ)
      
  
  !---------------------------------------------------------------    
  ! 温位
  !---------------------------------------------------------------
  xz_PotTempBasicZ = &
    & xz_TempBasicZ * (PressSfc / xz_PressBasicZ) ** (GasRDry / CpDry) 
  xz_PotTempBasicZ = xz_Boundary_Cyc_Sym(xz_PotTempBasicZ)

  
  !---------------------------------------------------------------    
  ! エクスナー関数
  !---------------------------------------------------------------
  xz_ExnerBasicZ = xz_TempBasicZ / xz_PotTempBasicZ    
  xz_ExnerBasicZ = xz_Boundary_Cyc_Sym(xz_ExnerBasicZ)

  
  !---------------------------------------------------------------    
  ! 密度
  !---------------------------------------------------------------
  xz_DensBasicZ = &
    & PressSfc * (xz_ExnerBasicZ ** (CvDry / GasRDry)) &
    &  / (GasRDry * xz_PotTempBasicZ / xz_EffMolWtBasicZ)
  xz_DensBasicZ = xz_Boundary_Cyc_Sym(xz_DensBasicZ)
  

  !---------------------------------------------------------------    
  ! 音速
  !---------------------------------------------------------------
  xz_VelSoundBasicZ = &
    & sqrt ( &
    &   CpDry * GasRDry * xz_ExnerBasicZ * xz_PotTempBasicZ &
    &   / (CvDry * xz_EffMolWtBasicZ) &
    & )
  xz_VelSoundBasicZ = xz_Boundary_Cyc_Sym(xz_VelSoundBasicZ)

 
  !----------------------------------------------------------
  ! BasicSet モジュールに値を設定  配列
  !----------------------------------------------------------
  call BasicSetArray_Init(                                  &
    & xz_PressBasicZ,    xz_ExnerBasicZ, xz_TempBasicZ,     &
    & xz_PotTempBasicZ,  xz_DensBasicZ,  xz_VelSoundBasicZ, &
    & xz_MixRtBasicZ, xz_EffMolWtBasicZ )

end subroutine BasicEnv
