!= Program ArareEnergy
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare-energy.f90,v 1.1 2009-08-24 04:07:06 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! 非静力学モデル deepconv/arare 2 次データ作成プログラム
!
!
program ArareEnergy

  !----- モジュール読み込み ------

  !-----   型宣言, 文字列処理   ----
  use dc_types,       only : STRING
  use dc_string,      only : StoA

  !-----   メッセージ出力   -----
  use dc_message,     only: MessageNotify

  !  コマンドライン引数解釈
  use argset,        only : argset_init  

  !-----    管理モジュール   -----
  !  化学量計算モジュール
  use chemcalc, only: chemcalc_init, xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH
  use chemdata, only: chemdata_init

  !  入出力ファイル名管理モジュール
  use fileset,       only : fileset_init

  !  デバッグ出力管理モジュール
  use debugset,      only : debugset_init

  !  時間管理モジュール
  use timeset,       only : timeset_init, TimeDisp, TImeInt

  !  格子点管理モジュール 
  use gridset,       only : gridset_init, DelX, DelZ, XMin, XMax, &
    &                       DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
    &                       RegXMax, RegZMax, RegXMin, RegZMin, s_Z

  !  基本場設定モジュール
  use basicset,      only : basicset_init,  &
    &                       MolWtDry, MolWtWet, CpDry, PressBasis, GasRDry,    &
    &                       SpcWetID, SpcWetMolfr, PressSfc, &
    &                       xza_MixRtBasicZ, xz_PotTempBasicZ, xz_ExnerBasicZ, &
    &                       xz_TempBasicZ, xz_PressBasicZ, Grav, xz_EffMolWtBasicZ, xz_DensBasicZ, xz_VelSoundBasicZ

  !  湿潤ルーチン設定モジュール
  use moistset,      only: moistset_init,   &
    &                      GasNum, CloudNum, RainNum, CondNum, RactNum, &
    &                      IdxG,   IdxC,   IdxR,                        &
    &                      IdxCG,  IdxCC,  IdxCR,                       &
    &                      IdxNH3, IdxH2S, IdxNH4SHc, IdxNH4SHr         

  use ChemCalc, only: xz_LatentHeat,    &!潜熱
    &                 ReactHeatNH4SH     !NH4SH の反応熱

  !  積算値管理モジュール
  use StorePotTemp,      only : StorePotTemp_init
  use StoreMixRt,        only : StoreMixRt_init
  use StoreBuoy,         only : StoreBuoy_init
  use StoreStab,         only : StoreStab_init

  !-----    下請けモジュール   -----
  !微分平均演算
  use average

  !  断熱上昇気塊の温度減率計算用モジュール
  use ECCM,         only : eccm_molfr
  
  ! 境界条件
  use boundary,      only : BoundaryXCyc_xz, BoundaryZSym_xz

  !  湿潤気塊の浮力計算用モジュール
  use MoistBuoyancy,only : MoistBuoy_Init, xr_BuoyMolWt, xr_BuoyDrag
  use DynFunc,       only : xr_Buoy

  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  real(8), allocatable :: xz_PotTemp(:,:)
  real(8), allocatable :: xz_PotTempAll(:,:)
  real(8), allocatable :: xz_Exner(:,:)
  real(8), allocatable :: pz_VelX(:,:)
  real(8), allocatable :: xr_VelZ(:,:)
  real(8), allocatable :: xza_MixRt(:,:,:)
  real(8), allocatable :: xza_MixRtAll(:,:,:)
  real(8), allocatable :: xza_MixRtAll2(:,:,:)
  real(8), allocatable :: xza_MixRtSat(:,:,:)
  real(8), allocatable :: xza_MixRtSatBasicZ(:,:,:)
  real(8)              :: AnalTime
  real(8)              :: Hum

  real(8), allocatable :: xz_MassDens(:,:)
  real(8), allocatable :: xz_KineticEnergy(:,:)
  real(8), allocatable :: xz_PotentialEnergy(:,:)
  real(8), allocatable :: xz_ElasticEnergyFO(:,:)
  real(8), allocatable :: xz_ElasticEnergySO(:,:)
  real(8)              :: MassTotal
  real(8)              :: KineticEnergyTotal
  real(8)              :: PotentialEnergyTotal
  real(8)              :: ElasticEnergyFOTotal
  real(8)              :: ElasticEnergySOTotal
!  real(8)              :: CloudMassTotal

  real(8), allocatable :: xz_TempAll(:,:)
  real(8), allocatable :: xz_PressAll(:,:)
  real(8), allocatable :: xz_EffMolWt(:,:)
  real(8), allocatable :: xza_MixRtDivMolWt(:,:,:)
  real(8), allocatable :: xza_LatentHeat(:,:,:)

  real(8), allocatable :: xz_MoistCond(:,:)
  real(8), allocatable :: xz_MoistSat(:,:)
  real(8), allocatable :: xz_MoistCondBasicZ(:,:)
  real(8), allocatable :: xz_MoistNH4SH(:,:)
  real(8), allocatable :: xz_MoistNH4SH2(:,:)
  real(8), allocatable :: xz_MoistNH4SH3(:,:)
  real(8), allocatable :: xz_DryENG(:,:)
  real(8), allocatable :: xz_DryEngBasicZ(:,:)
  real(8), allocatable :: xz_MoistENG(:,:)
  real(8), allocatable :: xz_SatMoistENG(:,:)
  real(8), allocatable :: xz_MoistENGBasicZ(:,:)

  real(8), allocatable :: xz_Z(:,:) 

  integer             :: t, s, i, k, tt, step
  character(100)       :: cfgfile

  !コマンドライン引数の解釈
  !  NAMELIST ファイル名の読み込み
  call argset_init(cfgfile)

  !物質特性の初期化
  call chemdata_init()

  !時刻に関する設定の初期化
  !  NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
  call timeset_init(cfgfile)

  !格子点情報の初期化
  !  NAMELIST から情報を得て, 格子点を計算する
  call gridset_init(cfgfile)
  
  !化学計算ルーチンの初期化
  call chemcalc_init()

  !基本場の情報の初期化
  !  NAMELIST から情報を得て, 基本場を設定する.
  call basicset_init(cfgfile)
  
  !I/O ファイル名の初期化
  !  NAMELIST ファイル名を指定し, deepconv/arare の
  !  出力ファイル名を NAMELIST から得る
  call fileset_init(cfgfile)
  
  !基本場の取得
  call AnalFile_BasicZ_Get( )

  !湿潤ルーチンの共有変数の初期化
  call moistset_init()

  !積算値を保管するためのモジュールの初期化
  !  NAMELIST から情報を得て, 基本場を設定する.
  call StorePotTemp_init( )
  call StoreMixRt_init( )
  call StoreBuoy_init( )  
  call StoreStab_init( )  

  !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
  call ArareAlloc

  !分子量に対する浮力計算ルーチンの初期化
  call MoistBuoy_Init()             
  
  !NetCDF ファイルの作成と基本場の取得  
!  call AnalFile_EqPotTemp_Open(  )
  call AnalFile_MoistENG_Open(  )
    
  !----------------------------------------------------------------------
  ! 解析値の計算開始
  !   時刻を進めながら, ファイルの値を得る
  !----------------------------------------------------------------------
!  do t = 1, int( TimeInt / TimeDisp ) + 1
  step = int( TimeInt / TimeDisp )
  STEPS: do while (step >= 600)
   write(*,*) step
   if (mod(step, 2) == 1) then 
       exit STEPS
    end if
    step = step / 2 
  end do STEPS
  write(*,*) "**** STEPS ***** ", int( TimeInt / TimeDisp ), " ==> ", step

  TIMESTEP: do tt = 0, int( TimeInt / TimeDisp) -1
     if ( mod(tt, int(TimeInt / TimeDisp / step)) /= 0 ) then 
        cycle TIMESTEP
     end if
     t = tt + 1
    
    !----------------------------------------------------------------
    ! ヒストリファイルを開き値を得る. 
    !----------------------------------------------------------------
    call AnalFile_Get(                                           &
      &       t,  AnalTime,                                      &
      &       xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt )
    write(*,*) t, AnalTime    
    
    !----------------------------------------------------------------
    ! 初期化
    !----------------------------------------------------------------
    xza_MixRtAll = xza_MixRt + xza_MixRtBasicZ

    !----------------------------------------------------------------
    ! 温度・圧力の計算
    !----------------------------------------------------------------
    xz_PressAll = &
      & PressBasis * ((xz_Exner + xz_ExnerBasicZ ) ** (CpDry / GasRDry))
    
    xz_PotTempAll  = xz_PotTemp + xz_PotTempBasicZ
    xz_TempAll = ( xz_Exner + xz_ExnerBasicZ ) * xz_PotTempAll

    !----------------------------------------------------------------
    ! 密度, 仮温位の計算
    !----------------------------------------------------------------
    xza_MixRtDivMolWt = 0.0d0
    do s = 1, GasNum
      xza_MixRtDivMolWt(:,:,IdxG(s)) = xza_MixRtAll(:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do

    xz_EffMolWt = &
      & (1.0d0 + sum(xza_MixRtAll,3)) &
      & / (MolWtDry * ((1.0d0 / MolWtDry) + sum(xza_MixRtDivMolWt,3))) 

    !----------------------------------------------------------------
    ! 飽和蒸気圧と平衡定数
    !----------------------------------------------------------------
    do s = 1, CondNum
      !飽和蒸気圧
      xza_MixRtSat(:,:,IdxCG(s)) =                        &
        & xz_SvapPress(SpcWetID(IdxCC(s)), xz_TempAll)    &
        &  * MolWtWet(IdxCC(s))                           &
        &  / (MolWtDry * xz_PressAll)                     

      xza_MixRtSatBasicZ(:,:,IdxCG(s)) =                  &
        & xz_SvapPress(SpcWetID(IdxCC(s)), xz_TempBasicZ) &
        &  * MolWtWet(IdxCC(s))                           &
        &  / (MolWtDry * xz_PressBasicZ)                  
    end do
    
    !---------------------------------------------------------------
    ! 湿潤静的エネルギー
    !---------------------------------------------------------------
    ! 鉛直座標を 2 次元配列化
    do k = DimZMin, DimZMax
      xz_Z(:,k) = s_Z(k)
    end do
    
    ! 乾燥静的エネルギー CpT+gz
    xz_DryEng = 0.0d0    
    xz_DryEngBasicZ = 0.0d0    
!    xz_DryEng       = CpDry * xz_TempAll + Grav / xz_EffMolWt * xz_Z
    xz_DryEng       = CpDry * xz_TempAll + Grav * xz_Z
    xz_DryEngBasicZ = CpDry * xz_TempBasicZ + Grav * xz_Z

    !NH4SH の生成反応. H2S は全て NH3 と反応して NH4SH になると考える. 
    ! H2S と反応した分だけ NH3 を減少させる. 
    xz_MoistNH4SH = ReactHeatNH4SH * xza_MixRtAll(:,:,IdxH2S) * MolWtWet(IdxNH4SHc) / MolWtWet(IdxH2S)
    xza_MixRtAll(:,:,IdxNH3) = xza_MixRtAll(:,:,IdxNH3) -  xza_MixRtAll(:,:,IdxH2S) &  
         &          * MolWtWet(IdxNH3) / MolWtWet(IdxH2S)  

    xz_MoistNH4SH2 = ReactHeatNH4SH * xza_MixRtSat(:,:,IdxH2S) * MolWtWet(IdxNH4SHc) / MolWtWet(IdxH2S)
    xza_MixRtSat(:,:,IdxNH3) = xza_MixRtSat(:,:,IdxNH3) -  xza_MixRtSat(:,:,IdxH2S) &  
         &          * MolWtWet(IdxNH3) / MolWtWet(IdxH2S)  

    xz_MoistNH4SH3 = ReactHeatNH4SH * xza_MixRtSatBasicZ(:,:,IdxH2S) * MolWtWet(IdxNH4SHc) / MolWtWet(IdxH2S)
    xza_MixRtSatBasicZ(:,:,IdxNH3) = xza_MixRtSatBasicZ(:,:,IdxNH3) - xza_MixRtSatBasicZ(:,:,IdxH2S) &  
         &          * MolWtWet(IdxNH3) / MolWtWet(IdxH2S)  

    ! H2O と NH3 の凝結の潜熱. 
    xz_MoistCond = 0.0d0
    xz_MoistSat = 0.0d0    
    xz_MoistCondBasicZ = 0.0d0    
    
    ! 凝結
    do s = 1, CondNum
      xz_MoistCond = xz_MoistCond     &
        &             + xz_LatentHeat( SpcWetID(IdxCC(s)), xz_TempAll )   &
        &               * xza_MixRtAll(:,:,IdxCG(s))  
    end do
    s = 1
!    do s = 1, CondNum
       xz_MoistSat = xz_MoistSat                                          &
        &             + xz_LatentHeat( SpcWetID(IdxCC(s)), xz_TempAll )   &
        &               * xza_MixRtSat(:,:,IdxCG(s))  
!    end do
    do s = 1, CondNum
      xz_MoistCondBasicZ = xz_MoistCondBasicZ                              &
        &             + xz_LatentHeat( SpcWetID(IdxCC(s)), xz_TempBasicZ ) &
        &               * xza_MixRtSatBasicZ(:,:,IdxCG(s))  
    end do    
    
    xz_MoistEng      = xz_DryEng       + xz_MoistCond       + xz_MoistNH4SH
!    xz_SatMoistEng   = xz_DryEng       + xz_MoistSat        + xz_MoistNH4SH2
    xz_SatMoistEng   = xz_DryEng       + xz_MoistSat
    xz_MoistEngBasicZ= xz_DryEngBasicZ + xz_MoistCondBasicZ + xz_MoistNH4SH3
  
    !----------------------------------------------------
    ! エネルギーに関するものたち
    !----------------------------------------------------
    ! 各格子点の質量密度(擾乱成分の寄与)
    xz_MassDens =                                                   &
      & (xz_ExnerBasicZ + xz_Exner)**( (CpDry - GasRDry)/GasRDry )  &
      & / (xz_PotTempBasicZ + xz_PotTemp)                           &
      & - (xz_ExnerBasicZ )**( (CpDry - GasRDry)/GasRDry )          &
      & / xz_PotTempBasicZ
    
    ! 領域全体の気相の質量
    MassTotal = (Xmax - Xmin) * PressSfc /Grav                 &
      & * ( xz_ExnerBasicZ(RegXMin,RegZMin)**(CpDry/GasRDry)          &
      &     - xz_ExnerBasicZ(RegXMin,RegZMax)**(CpDry/GasRDry) )      &
      & + PressSfc * DelX * DelZ / GasRDry                            &
      & * sum( xz_MassDens(RegXMin:RegXMax,RegZMin:RegZMax))
    
!    ! 各格子点の運動エネルギー(密度を和の密度で評価)
!     xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry      &
!      & * (xz_ExnerBasicZ + xz_Exner )**((CpDry - GasRDry)/GasRDry) &
!      & / (xz_PotTempBasicZ + xz_PotTemp)                           &
!      & * ( xz_avr_pz(pz_VelX)**2.0d0 + xz_avr_xr(xr_VelZ)**2.0d0 &
!      &     )

    ! 各格子点の運動エネルギー(密度を基本場の密度で評価)
    xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry      &
      & * (xz_ExnerBasicZ)**((CpDry - GasRDry)/GasRDry)               &
      & / (xz_PotTempBasicZ)                                          &
      & * ( xz_avr_pz(pz_VelX)**2.0d0 + xz_avr_xr(xr_VelZ)**2.0d0 &
      &     )
    
    ! 領域全体の運動エネルギー
    KineticEnergyTotal =                                             &
      & sum( xz_KineticEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
    
!    ! 領域全体の雲の質量
!    CloudMassTotal = DelX * DelZ                                     &
!      & * sum( xz_MixRt(RegXMin:RegXMax,RegZMin:RegZMax,1:SpcNum) )
    
    ! 各格子点の位置エネルギー
    xz_PotentialEnergy = - xz_DensBasicZ * xz_PotTemp * xz_Z       &
      & / xz_PotTempBasicZ
    
    ! 領域全体の位置エネルギー
    PotentialEnergyTotal = DelX * DelZ                               &
      & * sum( xz_PotentialEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
    
    ! 各格子点の弾性エネルギー(1 次量)
    xz_ElasticEnergyFO = xz_DensBasicZ *                             &
      & (CpDry * xz_PotTempBasicZ / xz_VelSoundBasicZ )**2.0d0        &
      & * xz_ExnerBasicZ * xz_Exner
    
    ! 領域全体の弾性エネルギー(1 次量)
    ElasticEnergyFOTotal = DelX * DelZ                               &
      & * sum( xz_ElasticEnergyFO(RegXMin:RegXMax,RegZMin:RegZMax) )
    
    ! 各格子点の弾性エネルギー(2 次量)
    xz_ElasticEnergySO = 0.5d0 * xz_DensBasicZ *                     &
      & (CpDry * xz_PotTempBasicZ * xz_Exner / xz_VelSoundBasicZ )**2.0d0
    
    ! 領域全体の弾性エネルギー(2 次量)
    ElasticEnergySOTotal = DelX * DelZ                               &
      & * sum( xz_ElasticEnergySO(RegXMin:RegXMax,RegZMin:RegZMax) )
    
    !----------------------------------------------------------------    
    ! ファイル出力
    !----------------------------------------------------------------
    call AnalFile_MoistENG_OutPut( AnalTime )
    
 end do TIMESTEP
  
  !----------------------------------------------------------------    
  ! ファイルを閉じる
  !----------------------------------------------------------------    
  call AnalFile_Close

    
contains
  
  subroutine AnalFile_Close( )
    !
    !解析ファイルのクローズ
    !
    use gt4_history, only: HistoryClose
    use fileset, only: gt_hist

    call HistoryClose
    
  end subroutine AnalFile_Close

  
  subroutine AnalFile_MoistENG_Open( )
    !
    !解析ファイルの定義
    !
    
    use gt4_history, only: HistoryCreate, HistoryPut, HistoryAddVariable
    use fileset, only: gt_hist, exptitle, expsrc, expinst, HistoryFilePrefix    
    use gridset, only: FileNX, FileNZ, s_X, s_Z, &
      &                FileXMin, FileXMax, FileZMin, FileZMax
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    character(150)  :: File

    !ファイル名の指定
    File = trim(HistoryFilePrefix)// '_Energy.nc'
    write(*,*) "******* OUTPUT FILE: ", File

    !-----------------------------------------------------------
    ! ヒストリー作成
    !-----------------------------------------------------------
    call HistoryCreate(                              &
      & file = File,                                 &
      & title = exptitle,                            &
      & source = expsrc,                             &
      & institution = expinst,                       &
      & dims=(/'x','z','t'/),                        &
      & dimsizes=(/FileNX, FileNZ, 0/),              &
      & longnames=(/'X-coordinate',                  &
      &             'Z-coordinate',                  &
      &             'Time        '/),                &
      & units=(/'m','m','s'/), origin=0.0,           &
      & interval=0.0 )
    
    !-----------------------------------------------------------  
    ! 軸の出力
    !-----------------------------------------------------------
    call HistoryPut('x', s_X( FileXMin: FileXMax ))
    call HistoryPut('z', s_Z( FileZMin: FileZMax ))

    !-----------------------------------------------------------  
    ! 解析用の変数の出力
    !-----------------------------------------------------------  
    call HistoryAddVariable(                                &
      & varname='Cp', dims=(/'x','z','t'/),           &
      & longname='Mean Specific Heat',                     &
      & units='J/K/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='CpBasicZ', dims=(/'x','z','t'/),           &
      & longname='Specific Heat (BasicZ)',                     &
      & units='J/K/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='MoistCond', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy (Lq)',                     &
      & units='J/kg', xtype='double' )
    
    call HistoryAddVariable(                                &
      & varname='MoistSat', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy (Lq_eq)',                     &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='MoistCondBasicZ', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy (Lq_sat, BasicZ)',                     &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='MoistNH4SH', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy (NH4SH)',                     &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='DrySEngBasicZ', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy (gz, BasicZ)',                     &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='MoistSEng', dims=(/'x','z','t'/),           &
      & longname='Moist Static Energy',                     &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='DrySEng', dims=(/'x','z','t'/),             &
      & longname='Dry Static Energy',                       &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                                &
      & varname='SatMoistSEng', dims=(/'x','z','t'/),        &
      & longname='Saturated Moist Static Energy',           &
      & units='J/kg', xtype='double' )


    call HistoryAddVariable(                                &
      & varname='MoistSEngBasicZ', dims=(/'x','z','t'/),  &
      & longname='Saturated Moist Static Energy (BasicZ)',  &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='MassDens', dims=(/'x','z','t'/),          &
      & longname='Mass Density',                           &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='MassTotal', dims=(/'t'/),                 &
      & longname='Mass Density',                           &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='KineticEnergy', dims=(/'x','z','t'/),     &
      & longname='Kinetic Energy',                         &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='KineticEnergyTotal', dims=(/'t'/),        &
      & longname='Total Kinetic Energy',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='PotentialEnergy', dims=(/'x','z','t'/),   &
      & longname='Potential Energy',                       &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='PotentialEnergyTotal', dims=(/'t'/),      &
      & longname='Total Potential Energy',                 &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='ElasticEnergyFO', dims=(/'x','z','t'/),   &
      & longname='Elastic Energy',                         &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='ElasticEnergyFOTotal', dims=(/'t'/),      &
      & longname='Total Elastic Energy',                   &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='ElasticEnergySO', dims=(/'x','z','t'/),   &
      & longname='Elastic Energy',                         &
      & units='J/kg', xtype='double' )

    call HistoryAddVariable(                               &
      & varname='ElasticEnergySOTotal', dims=(/'t'/),      &
      & longname='Total Elastic Energy',                   &
      & units='J/kg', xtype='double' )

  end subroutine AnalFile_MoistENG_Open


  subroutine AnalFile_MoistENG_OutPut( Time )

    use gridset, only: FileXMin, FileXMax, FileZMin, FileZMax
    use gt4_history, only: HistoryPut


    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in) :: Time    
 
    !----------------------------------------------------------------
    ! 値を出力
    !----------------------------------------------------------------
    call HistoryPut( 't', Time ) 
    
    call HistoryPut( 'Cp',          &
      &            CpDry / xz_EffMolWt(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'CpBasicZ',          &
      &            CpDry / xz_EffMolWtBasicZ(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'MoistCond',          &
      &            xz_MoistCond(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'MoistSat',          &
      &            xz_MoistSat(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'MoistCondBasicZ',          &
      &            xz_MoistCondBasicZ(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'MoistNH4SH',          &
      &            xz_MoistNH4SH(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'DrySEngBasicZ',          &
      &            xz_DryEngBasicZ(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'DrySEng',          &
      &            xz_DryENG(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'MoistSEng',          &
      &            xz_MoistENG(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'SatMoistSEng',          &
      &            xz_SatMoistENG(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'MoistSEngBasicZ',          &
      &            xz_MoistENGBasicZ(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'MassDens',          &
      &            xz_MassDens(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'MassTotal', MassTotal)

    call HistoryPut( 'KineticEnergy',          &
      &            xz_KineticEnergy(FileXMin:FileXMax, FileZMin:FileZMax))

    call HistoryPut( 'KineticEnergyTotal', KineticEnergyTotal )

    call HistoryPut( 'PotentialEnergy',          &
      &            xz_PotentialEnergy(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'PotentialEnergyTotal', PotentialEnergyTotal )

    call HistoryPut( 'ElasticEnergyFO',          &
      &            xz_ElasticEnergyFO(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'ElasticEnergyFOTotal', ElasticEnergyFOTotal)

    call HistoryPut( 'ElasticEnergySO',          &
      &            xz_ElasticEnergySO(FileXMin:FileXMax, FileZMin:FileZMax) )

    call HistoryPut( 'ElasticEnergySOTotal', ElasticEnergySOTotal )
    
  end subroutine AnalFile_MoistENG_OutPut


!!!------------------------------------------------------------------------------!!!
  subroutine AnalFile_Get(                                      &
    &       i, AnalTime,                                        &
    &       xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt )
 
    use dc_string
    use gt4_history,   only : HistoryGet
    use fileset,       only : HistoryFile
    use gridset,       only : DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
      &                       FileXMin, FileXMax, FileZMin, FileZMax
    use basicset,      only : SpcWetSymbol
    use boundary,      only : BoundaryXCyc_xz, BoundaryZSym_xz,   &
      &                       BoundaryXCyc_xza, BoundaryZSym_xza, &
      &                       BoundaryXCyc_pz, BoundaryZSym_pz,   &
      &                       BoundaryXCyc_xr, BoundaryZAntiSym_xr

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer, intent(in)  :: i    
    real(8), intent(out) :: AnalTime
    real(8), intent(out) :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(out) :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax,1:SpcNum)
    character(30)        :: name               !変数名
    character(10)        :: step
    
    step = 't=^' // adjustl(toChar(i))

    !-------------------------------------------------------------
    !Get a Value from netCDF File
    !-------------------------------------------------------------
    name = "t"
    call HistoryGet( HistoryFile(1), name, AnalTime, step )

    !-------------------------------------------------------------    
    ! Get a Value from netCDF File 
    !-------------------------------------------------------------
    name = "Exner"
    call HistoryGet( HistoryFile(1), name, xz_Exner(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xz( xz_Exner )
    call BoundaryZSym_xz( xz_Exner )
    
    name = "PotTemp"
    call HistoryGet( HistoryFile(2), name, xz_PotTemp(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xz( xz_PotTemp )
    call BoundaryZSym_xz( xz_PotTemp )

    name = "VelX"
    call HistoryGet( HistoryFile(3), name, pz_VelX(FileXMin:FileXMax, FileZMin:FileZMax), range=step )    
    call BoundaryXCyc_pz( pz_VelX )
    call BoundaryZSym_pz( pz_VelX )
    
    name = "VelZ"
    call HistoryGet( HistoryFile(4), name, xr_VelZ(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
    call BoundaryXCyc_xr( xr_VelZ )
    call BoundaryZAntiSym_xr( xr_VelZ )

!    name = "Km"
!    call HistoryGet( HistoryFile(5), name, xz_Km(FileXMin:FileXMax, FileZMin:FileZMax), range=step )
!    xz_Km = xz_BoundaryXCyc_xz( xz_Km )
!    xz_Km = xz_BoundaryZSym_xz( xz_Km )

    do s = 1, SpcNum
      name = trim(SpcWetSymbol(s))
      call HistoryGet( HistoryFile(8+s), name, xza_MixRt(FileXMin:FileXMax, FileZMin:FileZMax, s), range=step )
    end do

    call BoundaryXCyc_xza( xza_MixRt )
    call BoundaryZSym_xza( xza_MixRt )

  end subroutine AnalFile_Get



  subroutine AnalFile_BasicZ_Get( )

    use dc_string
    use gt4_history,   only: HistoryGet
    use fileset,       only: HistoryFile
    use basicset,      only: BasicSetArray_Init, SpcWetSymbol
    use gridset,       only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, &
      &                      FileXMin, FileXMax, FileZMin, FileZMax
    use boundary,      only: BoundaryXCyc_xz, BoundaryZSym_xz,     &
      &                      BoundaryXCyc_xza,BoundaryZSym_xza, &
      &                      BoundaryXCyc_pz, BoundaryZSym_pz,     &
      &                      BoundaryXCyc_xr, BoundaryZAntiSym_xr

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8) :: xz_DensBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_PotTempBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_ExnerBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xza_MixRtBZ(DimXMin:DimXMax, DimZMin:DimZMax,1:SpcNum)
    real(8) :: xz_EffMolWtBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_PressBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_TempBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) :: xz_VelSoundBZ(DimXMin:DimXMax, DimZMin:DimZMax)
    character(30)        :: name               !変数名
    

    !-------------------------------------------------------------
    ! 基本場の取得
    !-------------------------------------------------------------
    name = "TempBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_TempBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_TempBZ )
    call BoundaryZSym_xz( xz_TempBZ )

    name = "PressBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_PressBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_PressBZ )
    call BoundaryZSym_xz( xz_PressBZ )

    name = "ExnerBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_ExnerBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_ExnerBZ )
    call BoundaryZSym_xz( xz_ExnerBZ )

    name = "PotTempBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_PotTempBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_PotTempBZ )
    call BoundaryZSym_xz( xz_PotTempBZ )

    name = "DensBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_DensBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_DensBZ )
    call BoundaryZSym_xz( xz_DensBZ )

    name = "VelSoundBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_VelSoundBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_VelSoundBZ )
    call BoundaryZSym_xz( xz_VelSoundBZ )

    name = "EffMolWtBasicZ"
    call HistoryGet( HistoryFile(7), name, xz_EffMolWtBZ(FileXMin:FileXMax, FileZMin:FileZMax) )
    call BoundaryXCyc_xz( xz_EffMolWtBZ )
    call BoundaryZSym_xz( xz_EffMolWtBZ )

    do s = 1, SpcNum
      name = trim(SpcWetSymbol(s))//'BasicZ'
      call HistoryGet( HistoryFile(7), name, xza_MixRtBZ(FileXMin:FileXMax, FileZMin:FileZMax, s) )
    end do

    call BoundaryXCyc_xza( xza_MixRtBZ )
    call BoundaryZSym_xza( xza_MixRtBZ )

    !----------------------------------------------------------
    ! BasicSet モジュールに値を設定
    !----------------------------------------------------------
    call BasicSetArray_Init(                       &
      & xz_PressBZ,   xz_ExnerBZ,   xz_TempBZ,     &
      & xz_PotTempBZ, xz_DensBZ,    xz_VelSoundBZ, &
      & xza_MixRtBZ,  xz_EffMolWtBZ )

  end subroutine AnalFile_BasicZ_Get


  subroutine ArareAlloc

    use gridset,       only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum

    !基本場, 擾乱場の取得.
    allocate( &
      &       xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax ),             &
      &       xz_PotTempAll(DimXMin:DimXMax, DimZMin:DimZMax ),          &
      &       xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax ),               &
      &       pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax ),                &
      &       xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax ),                &
      &       xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),    &
      &       xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtAll2(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ), &
      &       xza_MixRtSatBasicZ(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum ),  &
      &       xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax ),            &
      &       xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      !
      &       xz_EffMolWt(DimXMin:DimXMax, DimZMin:DimZMax ),           &
      &       xza_MixRtDivMolWt(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),&
      &       xza_LatentHeat(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum),   &
      !
      &       xz_MoistCond(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistSat(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistCondBasicZ(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistNH4SH(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistNH4SH2(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistNH4SH3(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_DryEngBasicZ(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistENG(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_SatMoistENG(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_MoistENGBasicZ(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_DryENG(DimXMin:DimXMax,DimZMin:DimZMax),   &
      !
      &       xz_MassDens(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_KineticEnergy(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_PotentialEnergy(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_ElasticEnergyFO(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_ElasticEnergySO(DimXMin:DimXMax,DimZMin:DimZMax),   &
      &       xz_Z(DimXMin:DimXMax,DimZMin:DimZMax)   &
      &   )

    xza_MixRtSat = 0.0d0
    xza_MixRtDivMolWt = 0.0d0

    xz_MoistCond = 0.0d0
    xz_MoistSat = 0.0d0
    xz_MoistCondBasicZ = 0.0d0
    xz_MoistNH4SH = 0.0d0
    xz_DryEngBasicZ = 0.0d0
    xz_MoistENG = 0.0d0
    xz_SatMoistENG = 0.0d0
    xz_MoistENGBasicZ = 0.0d0
    xz_DryENG = 0.0d0

  end subroutine ArareAlloc

end program ArareEnergy
