!= Subroutine BasicEnv_3d
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: basicenv_3d.f90,v 1.3 2008-06-19 17:01:10 odakker Exp $
! Tag Name::  $Name: arare4-20100305 $
! 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
!

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

  !モジュール読み込み
  use dc_types, only : DP
  use dc_message, only: MessageNotify

  use gridset_3d,only: DimXMin,       &!配列の X 方向の下限
    &                 DimXMax,       &!配列の X 方向の上限
    &                 DimYMax,       &!配列の Y 方向の上限
    &                 DimYMin,       &!配列の Y 方向の下限
    &                 DimZMax,       &!配列の Z 方向の上限
    &                 DimZMin,       &!配列の Z 方向の下限
    &                 RegZMax,       &!配列の Z 方向の上限
    &                 RegZMin,       &!配列の Z 方向の下限
    &                 SpcNum,        &!凝縮成分の数
    &                 z_Z,           &!スカラー格子点での高度
    &                 r_dz           !Z 方向の格子点間隔
  use basicset_3d,only: BasicSetArray_Init, &!
    &                 PressBasis,    &!温位の基準圧力
    &                 GasRDry,       &!乾燥成分の定圧比熱
    &                 CpDry,         &!乾燥成分の定圧比熱
    &                 CvDry,         &!乾燥成分の定積比熱
    &                 MolWtDry,      &!乾燥成分の分子量
    &                 Grav,          &!重力加速度
    &                 SpcWetMolFr,   &!凝縮成分の初期モル比
    &                 MolWtWet,      &!凝縮成分の分子量
    &                 EnvType,       &!基本場の設定
    &                 Tropopause,    &!対流圏界面高度
    &                 GasRUniv,      &!
    &                 Humidity,      &!基本場の相対湿度
    &                 TempStrat,     &!成層圏の温度 [k]
!    &                 TempSfc,       &!表面温度 [k]
!    &                 PressSfc,      &!表面気圧 [Pa]
    &                 Dhight          !重み関数のパラメータ [m]
  use xyz_module, only: BoundaryXCyc_xyz, &
    &                   BoundaryYCyc_xyz, &
    &                   BoundaryZSym_xyz 

  use ECCM_3d,    only: ECCM_Dry,    &     
    &                   ECCM_Wet
  
  !暗黙の型宣言禁止
  implicit none

  !変数の定義
  real(DP) :: xyz_DensBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyz_PressBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyz_ExnerBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyz_TempBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyz_PotTempBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyz_VelSoundBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: xyza_MixRtBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)
  real(DP) :: xyz_EffMolWtBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
  real(DP) :: z_TempBasicZ(DimZMin:DimZMax)
  real(DP) :: z_PressBasicZ(DimZMin:DimZMax)
  real(DP) :: MolFrIni(SpcNum)
  real(DP) :: xyza_MolFr(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)
  real(DP) :: za_MolFr(DimZMin:DimZMax, SpcNum)
  real(DP) :: xyza_MixRtDivMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
  real(DP) :: z_DTempDZ(DimZMin:DimZMax)
  real(DP) :: z_MolWtMean(DimZMin:DimZMax)
  real(DP) :: DTempDZ
  real(DP) :: Weight
  integer  :: i, j, k, s

  !---------------------------------------------------------------
  ! 配列の初期化
  !---------------------------------------------------------------
  xyz_ExnerBasicZ    = 0.0d0
  xyz_PressBasicZ    = 0.0d0
  xyz_TempBasicZ     = 0.0d0
  xyz_PotTempBasicZ  = 0.0d0
  xyz_VelSoundBasicZ = 0.0d0
  xyza_MixRtBasicZ   = 0.0d0
  xyz_EffMolWtBasicZ = 0.0d0
  z_TempBasicZ      = 0.0d0
  z_PressBasicZ     = 0.0d0
  za_MolFr          = 0.0d0


  !---------------------------------------------------------------
  ! EnvType を元に, 温度, 圧力, 組成を決める
  !---------------------------------------------------------------
  MolFrIni = SpcWetMolFr(1:SpcNum) 

  !場合分け. Dry なら乾燥断熱的に, Wet なら湿潤断熱的な初期場を決める
  select case(EnvType)
  case("Dry")
    call ECCM_Dry( MolFrIni, Humidity, z_TempBasicZ, &
      &         z_PressBasicZ, z_MolWtMean, za_MolFr )
  case("Wet")
    call ECCM_Wet( MolFrIni, Humidity, z_TempBasicZ, &
      &         z_PressBasicZ, z_MolWtMean, za_MolFr )
  end select

  do k = RegZMin+1, DimZMax     
    z_DTempDZ(k) = (z_TempBasicZ(k) - z_TempBasicZ(k-1)) / r_dz(k-1)
  end do
  
  ! 対流圏界面より上の扱い
  !   圏界面より上は等温大気とする. 等温位大気から等温大気への遷移は 
  !   tanh を用いてなめらかにつなぐ
  do k = RegZMin+1, DimZMax
    
    !重みつけの関数を用意. tanh を用いる
    Weight = ( tanh( (z_Z(k) - tropopause ) / Dhight ) + 1.0d0 ) * 5.0d-1
    
    !仮値として温度を計算する. 圏界面より上では TempStrat の等温大気に近づける
    z_TempBasicZ(k) = z_TempBasicZ(k) * ( 1.0d0 - Weight ) + TempStrat * Weight
    
    !温度減率が断熱温度減率より小さくならないように
    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 * z_MolWtMean(k) / ( DTempDZ * GasRUniv ) ) )
  end do

  !確認のため出力
  call MessageNotify( "M", "BasicEnv", "Basic State Atmospheric Profiles." )
  do k = RegZMin, DimZMax-1     
    write(*,*) "temp", k, z_Z(k), z_TempBasicZ(k), z_PressBasicZ(k)
  end do

  ! 3 次元配列に格納
  do j = DimYMin, DimYMax
    do i = DimXMin, DimXMax
      xyz_TempBasicZ(i,j,:)  = z_TempBasicZ(:)  
      xyz_PressBasicZ(i,j,:) = z_PressBasicZ(:)
    end do
  end do

  !境界条件
  call BoundaryXCyc_xyz( xyz_TempBasicZ )
  call BoundaryYCyc_xyz( xyz_TempBasicZ )
  call BoundaryZSym_xyz( xyz_TempBasicZ )

  call BoundaryXCyc_xyz( xyz_PressBasicZ)
  call BoundaryYCyc_xyz( xyz_PressBasicZ )
  call BoundaryZSym_xyz( xyz_PressBasicZ )


!  case("IsoThermal")
!    z_TempBasicZ  = TempSfc
!
!    z_PressBasicZ = PressSfc * exp ( - Grav / z_TempBasicZ / GasRDry * z_Z)
!
!  end select

  !---------------------------------------------------------------
  ! 混合比
  !---------------------------------------------------------------
  !水平方向には一様
  do i = DimXMin, DimXMax      
    do j = DimYMin, DimYMax
      xyza_MolFr(i,j,:,:) = za_MolFr
    end do
  end do

  !気相のモル比を混合比に変換
  do s = 1, SpcNum
    xyza_MixRtBasicZ(:,:,:,s) = xyza_MolFr(:,:,:,s) * MolWtWet(s) / MolWtDry
  end do

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

  !---------------------------------------------------------------
  ! 分子量の効果
  !---------------------------------------------------------------
  do s = 1, SpcNum
    xyza_MixRtDivMolWt(:,:,:,s) = xyza_MixRtBasicZ(:,:,:,s) / MolWtWet(s)
  end do
  
  xyz_EffMolWtBasicZ = &
    & (1.0d0 + sum(xyza_MixRtBasicZ,4) ) &
    & / ( MolWtDry * ((1.0d0 / MolWtDry) + sum(xyza_MixRtDivMolWt,4)) )

  !境界条件  
  call BoundaryXCyc_xyz( xyz_EffMolWtBasicZ )
  call BoundaryYCyc_xyz( xyz_EffMolWtBasicZ )
  call BoundaryZSym_xyz( xyz_EffMolWtBasicZ )

  !---------------------------------------------------------------    
  ! 温位
  !---------------------------------------------------------------
  xyz_PotTempBasicZ = &
    & xyz_TempBasicZ * (PressBasis / xyz_PressBasicZ) ** (GasRDry / CpDry) 


  !境界条件  
  call BoundaryXCyc_xyz( xyz_PotTempBasicZ )
  call BoundaryYCyc_xyz( xyz_PotTempBasicZ )
  call BoundaryZSym_xyz( xyz_PotTempBasicZ )

  !---------------------------------------------------------------    
  ! エクスナー関数
  !---------------------------------------------------------------
  xyz_ExnerBasicZ = xyz_TempBasicZ / xyz_PotTempBasicZ    

  !境界条件
  call BoundaryXCyc_xyz( xyz_ExnerBasicZ )
  call BoundaryYCyc_xyz( xyz_ExnerBasicZ )
  call BoundaryZSym_xyz( xyz_ExnerBasicZ )
  
  !---------------------------------------------------------------    
  ! 密度
  !---------------------------------------------------------------
!  xyz_DensBasicZ = &
!    & PressBasis * (xyz_ExnerBasicZ ** (CvDry / GasRDry)) &
!    &  / (GasRDry * xyz_PotTempBasicZ )

  xyz_DensBasicZ = &
    & PressBasis * (xyz_ExnerBasicZ ** (CvDry / GasRDry)) &
    &  / (GasRDry * xyz_PotTempBasicZ / xyz_EffMolWtBasicZ)


  !境界条件
  call BoundaryXCyc_xyz( xyz_DensBasicZ )
  call BoundaryYCyc_xyz( xyz_DensBasicZ )
  call BoundaryZSym_xyz( xyz_DensBasicZ )

  !---------------------------------------------------------------    
  ! 音速
  !---------------------------------------------------------------
!  xyz_VelSoundBasicZ = &
!    & sqrt ( &
!    &   CpDry * GasRDry * xyz_ExnerBasicZ * xyz_PotTempBasicZ &
!    &   / CvDry &
!    & )

  xyz_VelSoundBasicZ = &
    & sqrt ( &
    &   CpDry * GasRDry * xyz_ExnerBasicZ * xyz_PotTempBasicZ &
    &   / (CvDry * xyz_EffMolWtBasicZ) &
    & )

  !境界条件
  call BoundaryXCyc_xyz( xyz_VelSoundBasicZ )
  call BoundaryYCyc_xyz( xyz_VelSoundBasicZ )
  call BoundaryZSym_xyz( xyz_VelSoundBasicZ )

  !----------------------------------------------------------
  ! BasicSet モジュールに値を設定
  !----------------------------------------------------------
!  call BasicSetArray_Init(                                  &
!    & xyz_PressBasicZ,    xyz_ExnerBasicZ, xyz_TempBasicZ,     &
!    & xyz_PotTempBasicZ,  xyz_DensBasicZ,  xyz_VelSoundBasicZ &
!    & )

  call BasicSetArray_Init(                                  &
    & xyz_PressBasicZ,    xyz_ExnerBasicZ, xyz_TempBasicZ,     &
    & xyz_PotTempBasicZ,  xyz_DensBasicZ,  xyz_VelSoundBasicZ, &
    & xyza_MixRtBasicZ, xyz_EffMolWtBasicZ )


end subroutine BasicEnv_3d
