!= Subroutine BasicEnv
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: basicenv_c2003.f90,v 1.1 2012-02-02 03:24:49 yamasita Exp $
! Tag Name::  $Name: arare4-20120511 $
! 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()
  !
  !デフォルトの基本場を設定するためのサブルーチン.
  !基本場を計算し, BasicSet モジュールの値を初期化する.
  !
  !コンパイルの順序の問題から, 基本場の値(hogeBasicZ な変数)を
  !計算する部分をBasicSet モジュールから切り離している.
  !ECCM 始め, BasicSet 自体に依存するが hogeBasicZ は use しない
  !外部サブルーチンを利用するためである.
  !

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

  use gridset,  only: DimXMin,       &!配列の X 方向の下限
    &                 RegXMin,       &!配列の X 方向の下限
    &                 DimXMax,       &!配列の X 方向の上限
    &                 RegXMax,       &!配列の X 方向の上限
    &                 DimZMin,       &!配列の Z 方向の下限
    &                 RegZMin,       &!配列の Z 方向の下限
    &                 DimZMax,       &!配列の Z 方向の上限
    &                 RegZMax,       &!配列の X 方向の上限
    &                 SpcNum,        &!凝縮成分の数
    &                 s_Z,           &!スカラー格子点での高度
    &                 DelZ            !Z 方向の格子点間隔
  use basicset, only: BasicSetArray_Init, &!
    &                 PressBasis,    &!温位の基準圧力
    &                 GasRDry,       &!乾燥成分の定圧比熱
    &                 CpDry,         &!乾燥成分の定圧比熱
    &                 CvDry,         &!乾燥成分の定積比熱
    &                 MolWtDry,      &!乾燥成分の分子量
    &                 Grav,          &!重力加速度
    &                 SpcWetMolFr,   &!凝縮成分の初期モル比
    &                 MolWtWet,      &!凝縮成分の分子量
    &                 EnvType,       &!基本場の設定
    &                 Tropopause,    &!対流圏界面高度
    &                 GasRUniv,      &!
    &                 Humidity,      &!基本場の相対湿度
    &                 TempStrat,     &!成層圏の温度 [k]
    &                 Dhight,        & !重み関数のパラメータ [m]
    &                 TempSfc,       &
    &                 PressSfc
  use Boundary, only: BoundaryXCyc_xz,  &!
    &                 BoundaryZSym_xz,  &!
    &                 BoundaryXCyc_xza, &!
    &                 BoundaryZSym_xza   !
  use ECCM,     only: ECCM_Dry,         &!
    &                 ECCM_Wet
  use cloudset, only: SatRtWetAdia
  use ChemData, only: ChemData_SVapPress_AntoineA, &
    &                 ChemData_SVapPress_AntoineB

  !暗黙の型宣言禁止
  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) :: xza_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  real(8) :: xz_EffMolWtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: z_TempBasicZ(DimZMin:DimZMax)
  real(8) :: z_PressBasicZ(DimZMin:DimZMax)
  real(8) :: MolFrIni(SpcNum)
  real(8) :: xza_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  real(8) :: za_MolFr(DimZMin:DimZMax, SpcNum)
  real(8) :: xza_MixRtDivMolWt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
  real(8) :: z_DTempDZ(DimZMin:DimZMax)
  real(8) :: z_MolWtMean(DimZMin:DimZMax)
  real(8) :: DTempDZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: Weight(DimXMin:DimXMax, DimZMin:DimZMax)
!  real(8) :: Weight1(DimXMin:DimXMax, DimZMin:DimZMax)
!  real(8) :: Weight2(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: LapseRateRatio ! 乾燥断熱減率に対する下層での温度減率の比
                            ! 乾燥断熱よりも小さく設定する.
  real(8) :: xz_TempAdia(DimXMin:DimXMax, DimZMin:DimZMax)
                                     ! 乾燥断熱線に沿った温度
  real(8) :: xz_TempSat(DimXMin:DimXMax, DimZMin:DimZMax)
                                     ! 凝結線に沿った温度
  real(8) :: xz_TempIso(DimXMin:DimXMax, DimZMin:DimZMax)
                                     ! 等温線に沿った温度
  real(8) :: xz_PressAdia(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_PressSat(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_PressIso(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_TempWork(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_PressWork(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8) :: xz_Z(DimXMin:DimXMax, DimZMin:DimZMax)
                                     ! 2D 座標
  real(8) :: Temp_0, Temp_1, Press_0, Press_1
                 ! 乾燥断熱線と湿潤断熱線とが交わる高度を反復法で求める
                 ! 際に用いる作業変数
  real(8) :: Work ! 作業変数
  real(8) :: LCL ! 下層の温度分布と湿潤断熱線が交わる高度
  real(8) :: LTP ! 湿潤断熱線と等温線が交わる高度
!  real(8) :: xz_LCL(DimXMin:DimXMax, DimZMin:DimZMax)
!             ! 下層の温度分布と湿潤断熱線が交わる高度
!  real(8) :: xz_LTP(DimXMin:DimXMax, DimZMin:DimZMax)
!             ! 下層の温度分布と湿潤断熱線が交わる高度
  integer :: i, k, s

  !---------------------------------------------------------------
  ! 配列の初期化
  !---------------------------------------------------------------
  xz_PressBasicZ    = 0.0d0
  xz_ExnerBasicZ    = 0.0d0
  xz_TempBasicZ     = 0.0d0
  xz_PotTempBasicZ  = 0.0d0
  xz_VelSoundBasicZ = 0.0d0
  xza_MixRtBasicZ   = 0.0d0
  xz_EffMolWtBasicZ = 0.0d0
!  z_TempBasicZ      = 0.0d0
  z_TempBasicZ      = TempSfc
  z_PressBasicZ     = 0.0d0
  za_MolFr          = 0.0d0

  ! 座標の初期化
    do k = DimZMin, DimZMax
       xz_Z(:,k) = s_Z(k)
    end do

  !---------------------------------------------------------------
  ! 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

! Colaprete et al.(2003) と同様, 飽和比一定であるような温度分布を与える.
! 但し上層は等温分布とする.

  if (TempSfc >= TempStrat ) then

       xz_TempIso  = TempStrat

       !--- 飽和比一定の温度分布と等温線が交わる高度(LTP)を計算

       LTP = - ChemData_SVapPress_AntoineB(12) * GasRDry / Grav &
         & * dlog( xz_TempIso(1,1) / ChemData_SVapPress_AntoineB(12) &
         & * ( ChemData_SVapPress_AntoineA(12) - dlog(PressSfc/SatRtWetAdia) &
         & ) )

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

       ! 飽和比一定の温度分布
       xz_TempSat = &
         & ChemData_SVapPress_AntoineB(12) / &
         & ( ChemData_SVapPress_AntoineA(12) - dlog(PressSfc/SatRtWetAdia) &
         & ) * &
         & exp( -Grav * xz_Z /(GasRDry*ChemData_SVapPress_AntoineB(12)) )

!        xz_TempSat = &
!             &  Temp_1 * exp(-Grav * (xz_Z - LCL) &
!             & /(GasRDry*ChemData_SVapPress_AntoineB(12)))

        Work = 3.0d3

!        Weight1 = ( tanh( (xz_Z - Work ) / Dhight ) + 1.0d0 ) * 5.0d-1
        Weight = ( tanh( (xz_Z - LTP ) / Dhight ) + 1.0d0 ) * 5.0d-1

        ! 温度勾配の不連続をなくす為に, tanh を用いてなめらかにつなぐ.

!        xz_TempWork = &
!             & xz_TempAdia * ( 1.0d0 - Weight1 ) &
!             & + xz_TempSat * Weight1
        xz_TempBasicZ = &
             & xz_TempSat * ( 1.0d0 - Weight ) + xz_TempIso * Weight
!        xz_TempBasicZ = &
!             & xz_TempWork * ( 1.0d0 - Weight2 ) + xz_TempIso * Weight2 &
!             & - 0.92d0

     else
        stop 'Tha value of TempSfc is invalid.'
     end if

  ! 圧力を静水圧平衡から計算
        do i = DimXMin, DimXMax
           DTempDZ(i,RegZMin) = &
                & (xz_TempBasicZ(i,RegZMin) - TempSfc) &
                & / ( 0.5d0 * DelZ)
           xz_PressBasicZ(i,RegZMin) =           &
                & Presssfc *  &
                & ( ( Tempsfc / xz_TempBasicZ(i,RegZMin) )  &
                &  ** (Grav / ( DTempDZ(i,RegZMin) * GasRDry ) ) )
        end do

        do k = RegZMin+1, DimZMax-1
           do i = DimXMin, DimXMax
              ! 局所的な温度減率
              DTempDZ(i,k) = (xz_TempBasicZ(i,k) - xz_TempBasicZ(i,k-1)) / DelZ
              xz_PressBasicZ(i,k) =                                   &
                   & xz_PressBasicZ(i,k-1) *                          &
                   & ( ( xz_TempBasicZ(i,k-1) / xz_TempBasicZ(i,k) )  &
                   &  ** (Grav / ( DTempDZ(i,k) * GasRDry ) ) )
           end do
        end do

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

  ! 2 次元配列に格納
!  do i = DimXMin, DimXMax
!    xz_TempBasicZ(i,:)  = z_TempBasicZ
!    xz_PressBasicZ(i,:) = z_PressBasicZ
!  end do

  !境界条件
  call BoundaryXCyc_xz( xz_TempBasicZ )
  call BoundaryZSym_xz( xz_TempBasicZ )
  call BoundaryXCyc_xz( xz_PressBasicZ )
  call BoundaryZSym_xz( xz_PressBasicZ )

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

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

!  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRtBasicZ <= 1.0d-20 )
!    xza_MixRtBasicZ = 1.0d-20
!  end where

  !境界条件
  call BoundaryXCyc_xza( xza_MixRtBasicZ )
  call BoundaryZSym_xza( xza_MixRtBasicZ )

  !---------------------------------------------------------------
  ! 分子量の効果
  !---------------------------------------------------------------
  do s = 1, SpcNum
    xza_MixRtDivMolWt(:,:,s) = xza_MixRtBasicZ(:,:,s) / MolWtWet(s)
  end do

  xz_EffMolWtBasicZ = &
    & (1.0d0 + sum(xza_MixRtBasicZ,3) ) &
    & / ( MolWtDry * ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3)) )

  !境界条件
  call BoundaryXCyc_xz( xz_EffMolWtBasicZ )
  call BoundaryZSym_xz( xz_EffMolWtBasicZ )

  !---------------------------------------------------------------
  ! 温位
  !---------------------------------------------------------------
  xz_PotTempBasicZ = &
    & xz_TempBasicZ * (PressBasis / xz_PressBasicZ) ** (GasRDry / CpDry)

  !境界条件
  call BoundaryXCyc_xz( xz_PotTempBasicZ )
  call BoundaryZSym_xz( xz_PotTempBasicZ )

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

  !境界条件
  call BoundaryXCyc_xz( xz_ExnerBasicZ )
  call BoundaryZSym_xz( xz_ExnerBasicZ )

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

  !境界条件
  call BoundaryXCyc_xz( xz_DensBasicZ )
  call BoundaryZSym_xz( xz_DensBasicZ )

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

  !境界条件
  call BoundaryXCyc_xz( xz_VelSoundBasicZ )
  call BoundaryZSym_xz( xz_VelSoundBasicZ )

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

end subroutine BasicEnv