!= Program Arare
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare_110725_falltest.f90,v 1.3 2011-10-26 16:55:48 yamasita Exp $
! Tag Name::  $Name: arare4-20111215 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! 非静力学モデル deepconv/arare. 
! 火星湿潤対流用. 
!
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 方程式系は準圧縮系.
!  * 温位の移流部分をひとまとめにして計算(08/12/02) 
!  * 位置エネルギー, 弾性エネルギーを出力(08/12/02)
!  * 雲密度の計算に乱流拡散, 数値粘性を導入(08/12/19)
!  * 圧力方程式に散逸加熱・放射加熱の効果を考慮(09/10/12)
!  * 質量収支チェック用プログラム. 
!    * 雲密度の計算には 
!      densitycloud_turb-masstest2.f90 (移流拡散だけ解く)
!      を用いること. 
!    * データ出力には historyfileio_mmconv2_masstest3.f90 を用いること. 
!    * 主成分凝結計算用穴埋めプログラムを導入. 
!  * 解析量の初期値を代入するよう書き換え. 
!  * 雲密度の 2 次元 tendency の計算手順を変更. 
!    * 短い時間ステップ内で和を取り, その平均量を保管
!    * 保管した値を出力時間間隔内で足しこむ
!    * 次の出力区間に移る前に値をクリア
!  * 雲密度の移流を長い時間ステップで評価
!  * やむなくゼロに引き戻した雲密度を出力するよう変更
!  * 雲密度の計算に時間フィルターを適用
!  * 凝結率の気相質量を全量で評価
!  * 雲粒落下を考慮したサーマル上昇実験用のテストプログラムファイル. 
!
!== Future Plans
!
!

program arare
  !
  !非静力学モデル deepconv/arare. 
  !

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

  !-----   型宣言, 文字列処理   ----
  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
  use chemdata, only: chemdata_init

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

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

  !  時間管理モジュール
  use timeset,       only : timeset_init, &
    &                       NstepLong, NstepShort, DelTimeLong, DelTimeShort, &
    &                       NstepDisp

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

  !  基本場設定モジュール
  use basicset,      only : basicset_init, &
    &                       xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, &
    &                       xz_VelSoundBasicZ, xz_ExnerBasicZ, &
    &                       CpDry, CvDry, GasRDry, Grav, PressSfc
 
  ! 平均操作モジュール(スカラー格子点, フラックス格子点間の変換)
  use average,       only : xz_avr_pz, xz_avr_xr, &
    &                       pz_avr_xz, xr_avr_xz

  ! 差分操作モジュール
  use differentiate_center4, only: xz_dx_pz, xz_dz_xr

  !  積算値管理モジュール
  use StorePotTemp,      only : StorePotTemp_init, StorePotTempClean, &
    &                           StorePotTempCond
  use StoreDensCloud,    only : StoreDensCloud_init, StoreDensCloudClean, &
    &                           StoreDensCloudAdv, StoreDensCloudTurb,  &
    &                           StoreDensCloudDiff, &
    &                           StoreDensCloudCond, StoreDensCloudFill, &
    &                           StoreDensCloudFillZero, &
    &                           StoreDensCloudMean, &
    &                           xz_DensCloudAdv,    &
    &                           xz_DensCloudTurb,   &
    &                           xz_DensCloudDiff,   &
    &                           xz_DensCloudCond,   &
    &                           xz_DensCloudFill,   &
    &                           xz_DensCloudFillZero,  &
    &                           xz_DensCloudFall,      &
    &                           xz_DensCloudAdv2,      &
    &                           xz_DensCloudTurb2,     &
    &                           xz_DensCloudDiff2,     &
    &                           xz_DensCloudCond2,     &
    &                           xz_DensCloudFill2,     &
    &                           xz_DensCloudFillZero2, &
    &                           xz_DensCloudFall2,     &
    &                           DensCloudAdvTendSum, DensCloudTurbTendSum, &
    &                           DensCloudDiffTendSum, DensCloudCondTendSum, &
    &                           DensCloudFillTendSum, DensCloudFillZeroTendSum, &
    &                           DensCloudFallTendSum
  use StoreMixRt,        only : StoreMixRt_init, StoreMixRtClean, &
    &                           StoreMixRtCond, StoreMixRtFill1, StoreMixRtFill2
  use StoreBuoy,         only : StoreBuoy_init, StoreBuoyClean
  use StoreStab,         only : StoreStab_init, StoreStabClean

  !  湿潤ルーチン設定モジュール
  use moistset,      only: moistset_init

  !-----    下請けモジュール   -----
  !  数値摩擦計算モジュール 
  use damping,       only : damping_init, & 
    &                        DampSponge_xz, DampSponge_pz, DampSponge_xr

  !  時間積分フィルターモジュール
  use timefilter,    only : AsselinFilter_xz, AsselinFilter_xr, &
    &                        AsselinFilter_pz, AsselinFilter_xza

  !  境界条件適用モジュール
  use boundary,    only : BoundaryXCyc_xz, BoundaryZSym_xz,     &
    &                      BoundaryXCyc_xza, BoundaryZSym_xza,   &
    &                      BoundaryXCyc_pz, BoundaryZSym_pz,     &
    &                      BoundaryXCyc_xr, BoundaryZAntiSym_xr

  !  CFL 条件確認モジュール
  use cflcheck,      only : CFLCheckTimeShort, &
    &                       CFLCheckTimeLongVelX, CFLCheckTimeLongVelZ

  !  負の湿潤量の補填計算モジュール
!  use fillnegative,  only : FillNegative_init, xza_FillNegative_xza
  use fillnegative,  only : xz_FillNegativeMMC_xz
  ! 主成分凝結対流用

  !-----    入出力モジュール   -----
  !  リスタートファイル入出力モジュール
  use RestartFileIO, only : ReStartFile_Open, ReStartFile_OutPut, &
    &                       ReStartFile_Close, ReStartFile_Get

  !  ヒストリファイル入出力モジュール
  use HistoryFileIO, only : HistoryFile_Open, HistoryFile_OutPut, &
    &                       HistoryFile_Close

  !-----       力学過程        -----
  !  力学過程計算用関数モジュール
  use DynFunc,       only : xz_AdvScalar, xz_AdvKm, xza_AdvScalar, pz_AdvVelX, &
    &                       xr_Buoy, xr_AdvVelZ, pz_GradPi
  
  !  力学過程陰解法計算用関数モジュール
!  use DynImpFunc,    only : xz_Exner_init, xz_Exner, xr_GradPi 
                                               ! 微量成分凝結用
  use DynImpFuncMarscond2,    only : xz_Exner_init, xz_Exner, xr_GradPi
                                               ! 主成分凝結用
  
  !-----       物理過程        -----
  !  数値拡散計算用モジュール
  use NumDiffusion,  only : NumDiffusion_Init, xz_NumDiffScalar, xz_NumDiffKm, &
    &                       xza_NumDiffScalar, pz_NumDiffVelX, xr_NumDiffVelZ

  !  乱流拡散計算用モジュール
  use Turbulence,   only : Turbulence_Init, &
    &                      xz_TurbScalar, xza_TurbScalar, pz_TurbVelX, &
    &                      xr_TurbVelZ  , xz_ShearKm    , xz_DispKm,   &
    &                      xz_DispHeat
  
  !  放射強制計算用モジュール
  use Radiation,    only : Radiation_init,  &
!    &                      xz_RadHeatBalance !加熱率と冷却率が各時刻で
!                                             !釣り合うよう調整
    &                      xz_RadHeatConst   !加熱率は一定

  
  !  地表フラックス計算用モジュール
!  use HeatFlux,     only : xz_HeatFluxBulk, xz_MixRtFluxBulk
!  use HeatFlux,     only : xz_HeatFluxDiff, xza_MixRtFluxDiff
  use HeatFlux_N1994,   only : xz_HeatFluxBulk, xz_MixRtFluxBulk

  !  湿潤飽和調節法計算用モジュール
  use MoistAdjust,  only : MoistAdjustSvapPress, MoistAdjustNH4SH

  !  雲物理パラメタリゼーション
  use WarmRainPrm,  only : WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, &
    &                      xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH,         &
    &                      xza_Cloud2Rain, xza_FallRain

  ! (2008/06/27 山下達也追加)
  ! 雲物理パラメータの初期化サブルーチン
  use cloudset, only : cloudset_init 


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

  !  安定度の計算
  use ECCM, only: ECCM_Stab


  !暗黙の型宣言禁止
  implicit none

  !内部変数
  character(80) :: cfgfile
  real(8), allocatable :: pz_VelXBl(:,:)
  real(8), allocatable :: pz_VelXNl(:,:)
  real(8), allocatable :: pz_VelXAl(:,:)
  real(8), allocatable :: pz_VelXNs(:,:)
  real(8), allocatable :: pz_VelXAs(:,:)
  real(8), allocatable :: xr_VelZBl(:,:)
  real(8), allocatable :: xr_VelZNl(:,:)
  real(8), allocatable :: xr_VelZAl(:,:)
  real(8), allocatable :: xr_VelZNs(:,:)
  real(8), allocatable :: xr_VelZAs(:,:)
  real(8), allocatable :: xz_ExnerBl(:,:)
  real(8), allocatable :: xz_ExnerNl(:,:)
  real(8), allocatable :: xz_ExnerAl(:,:)
  real(8), allocatable :: xz_ExnerNs(:,:)
  real(8), allocatable :: xz_ExnerAs(:,:)
  real(8), allocatable :: xz_ExnerSum(:,:)
  real(8), allocatable :: xz_PotTempWork(:,:)
  real(8), allocatable :: xz_PotTempBl(:,:)
  real(8), allocatable :: xz_PotTempNl(:,:)
  real(8), allocatable :: xz_PotTempAl(:,:)
  real(8), allocatable :: xz_PotTempNs(:,:)
  real(8), allocatable :: xz_PotTempAs(:,:)
  real(8), allocatable :: xz_PotTempSum(:,:)
  real(8), allocatable :: xz_TempSum(:,:)
  real(8), allocatable :: xz_KmBl(:,:)
  real(8), allocatable :: xz_KmNl(:,:)
  real(8), allocatable :: xz_KmAl(:,:)
  real(8), allocatable :: xz_KhBl(:,:)
  real(8), allocatable :: xz_KhNl(:,:)
  real(8), allocatable :: xz_KhAl(:,:)
  real(8), allocatable :: xza_MixRtWork(:,:,:)
  real(8), allocatable :: xza_MixRtBl(:,:,:)
  real(8), allocatable :: xza_MixRtNl(:,:,:)
  real(8), allocatable :: xza_MixRtAl(:,:,:)

  real(8), allocatable :: pz_AccelVelXNl(:,:)
  real(8), allocatable :: xr_AccelVelZNl(:,:)
  real(8), allocatable :: xza_DelMixRt(:,:,:)
  real(8), allocatable :: xz_FluxDensCloud(:,:)
  real(8), allocatable :: xz_TurbDensCloud(:,:)
  real(8), allocatable :: xz_NumDensCloud(:,:)

! 長い時間ステップで評価すべき温位の項を格納する配列
! 2008/06/20 山下達也 追加
  real(8), allocatable :: xz_TendDensCloudNl(:,:)
  real(8), allocatable :: xz_TendPotTempNl(:,:)

                         

  real(8)              :: Time
  real(8)              :: ReStartTime(2)
  real(8), allocatable :: DelTimeLFrog(:)
  real(8)              :: DelTimeEular
  integer, allocatable :: NStepEular(:)
  integer              :: NStepLFrog

  integer              :: t, tau, t1, t2, k

! 以下の変数は山下が追加(2008/05/07)

! 過飽和度
  real(8), allocatable :: xz_SatRatioBl(:,:)
  real(8), allocatable :: xz_SatRatioNl(:,:)
  real(8), allocatable :: xz_SatRatioAl(:,:)
  real(8), allocatable :: xz_SatRatioNs(:,:)
  real(8), allocatable :: xz_SatRatioAs(:,:)
! 雲密度
  real(8), allocatable :: xz_DensCloudBl(:,:)
  real(8), allocatable :: xz_DensCloudNl(:,:)
  real(8), allocatable :: xz_DensCloudAl(:,:)
  real(8), allocatable :: xz_DensCloudNs(:,:)
  real(8), allocatable :: xz_DensCloudAs(:,:)

  real(8), allocatable :: xz_DensCloudWorkAs(:,:) ! 作業変数
  real(8)              :: DensCloudMin ! 負の雲密度穴埋めの為の作業変数

! 凝結量
  real(8), allocatable :: xz_MassCondNs(:,:)
  real(8), allocatable :: xz_MassCondNl(:,:)
! 凝結熱による温位変化率
  real(8), allocatable :: xz_QCond(:,:)
! 放射による温位変化率
  real(8), allocatable :: xz_QRad(:,:)
! 熱散逸による温位変化率
  real(8), allocatable :: xz_QDis(:,:)
! 作業変数
  real(8), allocatable :: worknum(:,:)
  real(8)              :: workcloud
! 各格子点の質量密度(擾乱成分の寄与)
  real(8), allocatable :: xz_MassDens(:,:)
! 領域全体の気相の質量
  real(8)              :: MassTotal
! 温位 2 次移流項の全質量の時間変化に対する寄与
  real(8), allocatable :: xz_MassTend(:,:)
  real(8)              :: MassTend
! 領域全体の運動エネルギー
  real(8)              :: KineticEnergyTotal
! 各格子点の運動エネルギー
  real(8), allocatable :: xz_KineticEnergy(:,:)
! 領域全体の雲の質量
  real(8)              :: CloudMassTotal
! 領域全体の雲の質量(fillnegative の前, write 文でチェック)
  real(8)              :: CloudMassTotal01
! 領域全体の雲の質量(fillnegative の後, write 文でチェック)
  real(8)              :: CloudMassTotal02
! 各格子点のポテンシャルエネルギー
  real(8), allocatable :: xz_PotentialEnergy(:,:)
! 領域全体のポテンシャルエネルギー
  real(8)              :: PotentialEnergyTotal
! 各格子点の弾性エネルギー(1 次量)
  real(8), allocatable :: xz_ElasticEnergyFO(:,:)
! 領域全体の弾性エネルギー(1 次量)
  real(8)              :: ElasticEnergyFOTotal
! 各格子点の弾性エネルギー(2 次量)
  real(8), allocatable :: xz_ElasticEnergySO(:,:)
! 領域全体の弾性エネルギー(2 次量)
  real(8)              :: ElasticEnergySOTotal
! z 座標(2 次元配列)
  real(8), allocatable :: xz_Z(:,:)


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

  !デバッグ設定
  call debugset_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 cloudset_init(cfgfile)

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

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

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

  !予報変数の初期化
  !  ReStartFile が設定されている場合にはファイルを読み込み, 
  !  設定されていない場合には, デフォルトの基本場と擾乱場を作る. 

  write(*,*) InitFile

  if (trim(InitFile) /= '') then    

    !基本場, 擾乱場の初期値を netCDF ファイルから取得する.
    call ReStartFile_Get(            &
      & ReStartTime,                 &
      & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
      & xza_MixRtBl,  xz_KmBl,    xz_KhBl,              &
      & xz_DensCloudBl, xz_SatRatioBl,                  &
      & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
      & xza_MixRtNl,  xz_KmNl,    xz_KhNl,              &
      & xz_DensCloudNl, xz_SatRatioNl                    )

  else
    !デフォルト設定の基本場, 擾乱場を作成する. 
    call BasicEnv()
    call DisturbEnv(cfgfile,                            &
      & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
      & xza_MixRtBl,  xz_KmBl,    xz_KhBl,              &
      & xz_DensCloudBl, xz_SatRatioBl                     )
    
    ! 1 ループ目だけは, bl と nl の値を同じにしておく. 
    ! 1 ステップ目はオイラー法で解く必要があるが, 
    ! 1 ステップ目とそれ以外のステップを別々にコーディングしたくない為
    xz_PotTempNl = xz_PotTempBl
    xz_ExnerNl   = xz_ExnerBl
    pz_VelXNl    = pz_VelXBl
    xr_VelZNl    = xr_VelZBl
    xza_MixRtNl  = xza_MixRtBl
    xz_KmNl      = xz_KmBl
    xz_KhNl      = xz_KhBl
    xz_DensCloudNl = xz_DensCloudBl
    xz_SatRatioNl = xz_SatRatioBl
    xz_PotTempSum = xz_PotTempBl + xz_PotTempBasicZ
    xz_ExnerSum = xz_ExnerBl + xz_ExnerBasicZ
    xz_TempSum = xz_PotTempSum * xz_ExnerSum
    xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry       &
      & * (xz_ExnerBasicZ)**( CvDry/GasRDry )                         &
      & / (xz_PotTempBasicZ)                                          &
      & * ( xz_avr_pz(pz_VelXNl)**2.0d0 + xz_avr_xr(xr_VelZNl)**2.0d0 &
      &     )
    KineticEnergyTotal =                                              &
      & sum( xz_KineticEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
    xz_MassDens =                                                     &
      & (xz_ExnerBasicZ + xz_ExnerNl)**( CvDry/GasRDry )              &
      & / (xz_PotTempBasicZ + xz_PotTempNl)                           &
      & - (xz_ExnerBasicZ )**( CvDry/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))
    CloudMassTotal = DelX * DelZ                                     &
      & * sum( xz_DensCloudNl(RegXMin:RegXMax,RegZMin:RegZMax) )
    do k = RegZMin, RegZMax
       xz_Z(:,k) = s_Z(k)
    end do
    xz_PotentialEnergy = - xz_DensBasicZ * xz_PotTempNl * xz_Z       &
      & / xz_PotTempBasicZ
    PotentialEnergyTotal = DelX * DelZ                               &
      & * sum( xz_PotentialEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
    xz_MassTend = - xz_DensBasicZ * xz_TendPotTempNl / xz_PotTempBasicZ
    MassTend = DelX * DelZ                                     &
         & * sum( xz_MassTend(RegXMin:RegXMax,RegZMin:RegZMax) )

  end if

  !----------------------------------------------------------------------
  ! パッケージ型モジュールの初期化
  !   デフォルトの値から変更する必要のあるルーチンのみ初期化
  !----------------------------------------------------------------------
  call Damping_Init( cfgfile )      !波の減衰係数の初期化
  call NumDiffusion_Init()          !数値拡散項の初期化
  call Turbulence_Init()            !乱流計算の初期化
  call WarmRainPrm_Init( cfgfile )  !暖かい雨のパラメタリゼーションの初期化
!  call FillNegative_Init( xza_MixRtBasicZ, xz_DensBasicZ) 
!                                    !移流による負の量の処理
  call Radiation_Init( cfgfile )    !放射強制の初期化
  call MoistBuoy_Init()             !分子量に対する浮力計算ルーチンの初期化
  call xz_Exner_Init()              !陰解法の初期化  
       

  !----------------------------------------------------------------------
  ! 時刻とループ回数の初期化 
  !----------------------------------------------------------------------
  NstepLFrog   = NstepLong 
  NstepEular   = NstepShort 
  DelTimeLFrog = DelTimeLong * 2.0d0 
  DelTimeEular = DelTimeShort
  if ( trim(InitFile) /= '') then    
    Time = ReStartTime(2)               !リスタート開始時刻
  else
    Time = 0.0d0                          !計算開始時刻
    NstepEular(1)   = NstepShort /2 ! 1 ループ目だけ
    DelTimeLFrog(1) = DelTimeLong         ! 1 ループ目だけ
  end if

 
  !----------------------------------------------------------------
  ! ファイル入出力
  !----------------------------------------------------------------
  !ファイルオープン
  call HistoryFile_Open( )

  if ( trim(InitFile) /= '') then    
    call HistoryFile_Output(  &
      &  ReStartTime(2),      &
      &  xz_PotTempNl,        &
      &  xz_PotTempSum,       &
      &  xz_TempSum,          &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xza_MixRtNl,         &
      &  xz_KmNl,             &
      &  xz_KhNl,             &
      &  xz_DensCloudNl,      &
      &  xz_SatRatioNl,       &
      &  MassTotal,           &
      &  MassTend,            &
      &  KineticEnergyTotal,  &
      &  CloudMassTotal,      &
      &  PotentialEnergyTotal, &
      &  ElasticEnergyFOTotal, &
      &  ElasticEnergySOTotal   )
  end if

  !時刻 0 の場合には出力
  if ( Time == 0 ) then
    call HistoryFile_Output(  &
      &  Time,                &
      &  xz_PotTempNl,        &
      &  xz_PotTempSum,       &
      &  xz_TempSum,          &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xza_MixRtNl,         &
      &  xz_KmNl,             &
      &  xz_KhNl,             &
      &  xz_DensCloudNl,      &
      &  xz_SatRatioNl,       &
      &  MassTotal,           &
      &  MassTend,            &
      &  KineticEnergyTotal,  &
      &  CloudMassTotal,      &
      &  PotentialEnergyTotal, &
      &  ElasticEnergyFOTotal, &
      &  ElasticEnergySOTotal   )
  end if
  
  !----------------------------------------------------------------------
  ! 設定のチェック
  !----------------------------------------------------------------------
  !CFL 条件のチェック
  call CFLCheckTimeShort( xz_VelSoundBasicZ )

  !----------------------------------------------------------------------
  ! 数値積分
  !----------------------------------------------------------------------
  call MessageNotify( "M", "main", "Time integration" ) 
 
  do t1 = 1, NstepLFrog / NstepDisp
    do t2 = 1, NstepDisp

      t = (t1 - 1) * NstepDisp + t2
      
      !時刻の設定
      Time = Time + DelTimeLong

      !----------------------------------------------------------------
      ! 渦粘性係数, 渦拡散係数を求める.
      !----------------------------------------------------------------
      xz_KmAl =                                                       &
        & xz_KmBl                                                     &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)               &
        &     + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) &
        &     + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)             &   
        &     + xz_NumDiffKm(xz_KmBl)                                 &
        &     + xz_DispKm(xz_KmBl)                                    &
        &    )

      !値の上限下限の設定
      !  * 値は正になることを保証する
      !  * 値の上限は 800 とする. 中島健介(1994, 学位論文)参照
      xz_KmAl = max( 0.0d0, min( xz_KmAl, 800.0d0 ) )
!     xz_KmAl = max( 0.0d0, min( xz_KmAl, 300.0d0 ) )

      !境界条件
      call BoundaryXCyc_xz( xz_KmAl )
      call BoundaryZSym_xz( xz_KmAl )
      
      !渦拡散定数を決める
      xz_KhAl = 3.0d0 * xz_KmAl
      
      !-------------------------------------------------------------
      ! 長い時間ステップでの速度の移流, 拡散, 数値粘性, 浮力
      !                     温位の移流, 拡散, 数値粘性
      !                     雲密度の移流, 拡散, 数値粘性
      !-------------------------------------------------------------

      ! 雲密度のフラックス項の計算. 4 次精度中心差分を用いて計算

      xz_FluxDensCloud  =                                             &
        & - xz_dx_pz(pz_VelXNl * pz_avr_xz(xz_DensCloudNl))           &
        & - xz_dz_xr(xr_VelZNl * xr_avr_xz(xz_DensCloudNl))

      call StoreDensCloudAdv( xz_FluxDensCloud )

      ! 雲密度の乱流拡散項の計算

      xz_TurbDensCloud = xz_TurbScalar(xz_DensCloudBl, xz_KhBl)

      call StoreDensCloudTurb( xz_TurbDensCloud )

      ! 雲密度の数値粘性項の計算

      xz_NumDensCloud = xz_NumDiffScalar(xz_DensCloudBl)

      call StoreDensCloudDiff( xz_NumDensCloud )

      xz_TendDensCloudNl =                                            &
        &     xz_FluxDensCloud                                        &
        &     + xz_TurbDensCloud                                      &
        &     + xz_NumDensCloud


      xz_PotTempSum = xz_PotTempNl + xz_PotTempBasicZ

      xz_TendPotTempNl =                                              &
        &     + xz_AdvScalar( xz_PotTempSum,    pz_VelXNl, xr_VelZNl) &
        &     + xz_TurbScalar(xz_PotTempBl,     xz_KhBl)              &
        &     + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)              &
        &     + xz_NumDiffScalar(xz_PotTempBl)                        

      pz_AccelVelXNl =                                    &
        & + pz_AdvVelX(pz_VelXNl, xr_VelZNl)              &
        & + pz_TurbVelX(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
        & + pz_NumDiffVelX(pz_VelXBl)
      
      xr_AccelVelZNl =                                    &
        & + xr_AdvVelZ(xr_VelZNl, pz_VelXNl)              &
        & + xr_Buoy(xz_PotTempNl)                         &    
!        & + xr_BuoyMolWt(xza_MixRtNl)                     &
!        & + xr_BuoyDrag(xza_MixRtNl)                      &
        & + xr_TurbVelZ(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
        & + xr_NumDiffVelZ(xr_VelZBl)                        

      !-------------------------------------------------------------
      ! 短い時間ステップの初期値作成
      !-------------------------------------------------------------
      pz_VelXNs  = pz_VelXBl
      xr_VelZNs  = xr_VelZBl
      xz_ExnerNs = xz_ExnerBl
      xz_PotTempNs = xz_PotTempBl
      xz_DensCloudNs = xz_DensCloudBl
      xz_SatRatioNs = xz_SatRatioBl      


      !-------------------------------------------------------------
      ! 短い時間ステップの積分. オイラー法を利用
      !-------------------------------------------------------------
      Eular: do tau = 1, NstepEular(t)

      !-------------------------------------------------------------
      ! 2008/04/14 山下追加
      ! 雲密度の計算, 凝結熱の計算, 凝結量の計算
      !-------------------------------------------------------------

      !=== 雲密度の移流・拡散・落下計算
      call DensityCloud(  &
        & xz_DensCloudNs,                        & !(in) 現在の時間の雲密度
        & DelTimeShort,                          & !(in) 時間間隔
        & xz_DensCloudWorkAs )                     !(out)雲密度(負の密度調整前)

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

      xz_DensCloudAs = xz_FillNegativeMMC_xz( xz_DensCloudWorkAs ) 

      ! 埋めた/削った量を保管
      call StoreDensCloudFill( (xz_DensCloudAs - xz_DensCloudWorkAs) / DelTimeShort )

      xz_DensCloudWorkAs = xz_DensCloudAs
      xz_DensCloudAs = max( 0.0d0, xz_DensCloudAs )

      ! やむなくゼロに戻した量を保管
      call StoreDensCloudFillZero( (xz_DensCloudAs - xz_DensCloudWorkAs) / DelTimeShort )

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

      xz_MassCondNs = 0.0d0
     call StoreDensCloudCond( xz_MassCondNs )

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

      ! (2005/12/07 北守太一)
      !   * 凝結により発生した潜熱を計算
      !=== 凝結熱の計算
      call LatentHeat(  & 
        & xz_MassCondNs,       & !(in) 凝結量
        & xz_ExnerNs,          & !(in) エクスナー関数擾乱成分
        & xz_PotTempNs,        & !(in) 温位擾乱成分
        & xz_QCond)              !(out)凝結熱による温位変化率

      !
      ! 温位の移流計算. 
      !(2008/06/20, 山下達也 : 
      ! 主成分凝結を議論するため, 短い時間ステップのループ内
      ! で計算するよう書き換えた. )
      !


      xz_PotTempAs =                                                  &
        &   xz_PotTempNs                                              &
        & + DelTimeEular                                              &
        &   * (                                                       &
!        &     + xz_AdvScalar( xz_PotTempNl,     pz_VelXNl, xr_VelZNl) &
!        &     + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl) &
!        &     + xz_TurbScalar(xz_PotTempBl,     xz_KhBl)              &
!        &     + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)              &
!        &     + xz_NumDiffScalar(xz_PotTempBl)                        &
        &     + xz_TendPotTempNl                                      &
        &     + xz_DispHeat( xz_KmBl )                                &
!        &     + xz_RadHeatConst( xz_ExnerBasicZ )                      &
!        &    + xz_RadHeatBalance( xz_ExnerBasicZ, xz_PotTempBasicZ, xz_ExnerNs, xz_PotTempNs )         &
!        &     + xz_HeatFluxDiff( xz_PotTempNl )                       &
        &     + xz_Qcond                                              &
!      &     + xz_HeatFluxBulk( xz_PotTempNl, pz_VelXNl )             &
!      &     + xz_HeatFluxBulk( xz_PotTempNl  )                      &
!!      &     + xz_NewtonCool( xz_PotTempBl )                       &
        &      )

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



        
        !-------------------------------------------------------------
        ! 速度 u の計算
        !-------------------------------------------------------------
        pz_VelXAs = 0.0d0
        !境界条件
        call BoundaryXCyc_pz( pz_VelXAs )
        call BoundaryZSym_pz( pz_VelXAs )
        
        !-------------------------------------------------------------
        ! エクスナー関数の計算
        !-------------------------------------------------------------

        xz_QRad = 0.0d0
        xz_QDis = xz_DispHeat( xz_KmBl )

        xz_ExnerAs = 0.0d0
        
        !境界条件
        call BoundaryXCyc_xz( xz_ExnerAs )
        call BoundaryZSym_xz( xz_ExnerAs )
        
        !-------------------------------------------------------------
        ! 速度 w の計算
        !-------------------------------------------------------------
        xr_VelZAs = 0.0d0
        !境界条件
        call BoundaryXCyc_xr( xr_VelZAs )
        call BoundaryZAntiSym_xr( xr_VelZAs )

      ! (2008/06/16 山下達也)
      !=== 過飽和度の計算
      call SaturationRatio( &
        & xz_ExnerAs,       & !(in)
        & xz_PotTempAs,     & !(in)
        & xz_SatRatioAs)     !(out)

        
        !-------------------------------------------------------------
        ! 短い時間ステップのループを回すための処置
        !-------------------------------------------------------------
        xz_ExnerNs  = xz_ExnerAs
        xz_PotTempNs = xz_PotTempAs
        pz_VelXNs   = pz_VelXAs
        xr_VelZNs   = xr_VelZAs
        xz_DensCloudNs = xz_DensCloudAs
        xz_SatRatioNs = xz_SatRatioAs

      end do Eular

      ! 短いループでの tendency を平均
      call StoreDensCloudMean( ) 

      ! 平均された tendency の足しこみ
      call DensCloudAdvTendSum( xz_DensCloudAdv )
      call DensCloudTurbTendSum( xz_DensCloudTurb )
      call DensCloudDiffTendSum( xz_DensCloudDiff )
      call DensCloudCondTendSum( xz_DensCloudCond )
      call DensCloudFillTendSum( xz_DensCloudFill )
      call DensCloudFillZeroTendSum( xz_DensCloudFillZero )
      call DensCloudFallTendSum( xz_DensCloudFall )

      !----------------------------------------------------------------
      ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
      !----------------------------------------------------------------
      xz_ExnerAl  = xz_ExnerAs
      xz_PotTempAl = xz_PotTempAs
      pz_VelXAl   = pz_VelXAs
      xr_VelZAl   = xr_VelZAs
      xz_DensCloudAl = xz_DensCloudAs
      xz_SatRatioAl = xz_SatRatioAs


      !-------------------------------------------------------------
      ! アッセリンの時間フィルタ.  Nl の値をフィルタリング
      ! 1 ステップ目は Euler 法で計算するため, フィルタをかけない. 
      ! (2008/08/25 山下達也)
      !-------------------------------------------------------------

  if ( t /= 1 ) then
      call AsselinFilter_xz(xz_ExnerAl, xz_ExnerNl, xz_ExnerBl)
      call AsselinFilter_pz(pz_VelXAl, pz_VelXNl, pz_VelXBl )
      call AsselinFilter_xr(xr_VelZAl, xr_VelZNl, xr_VelZBl )
      call AsselinFilter_xz(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl)
      call AsselinFilter_xz(xz_KmAl, xz_KmNl, xz_KmBl)
      call AsselinFilter_xz(xz_DensCloudAl, xz_DensCloudNl, xz_DensCloudBl)
!      call AsselinFilter_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl)

  else
      write(*,*) "OK"
  end if
    
      !-------------------------------------------------------------
      ! スポンジ層
      !-------------------------------------------------------------
      call DampSponge_pz( pz_VelXAl, pz_VelXBl, DelTimeLFrog(t) )
      call DampSponge_xr( xr_VelZAl, xr_VelZBl, DelTimeLFrog(t) )
      call DampSponge_xz( xz_PotTempAl,  xz_PotTempBl, DelTimeLFrog(t) )
      call DampSponge_xz( xz_ExnerAl,    xz_ExnerBl,   DelTimeLFrog(t) )
      call DampSponge_xz( xz_KmAl,       xz_KmBl,      DelTimeLFrog(t) )

      !--------------------------------------------------------------
      ! 解析値
      !---------------------------------------------------------------
      call ECCM_Stab( xz_PotTempAl, xz_ExnerAl, xza_MixRtAl )
      
      !----------------------------------------------------------------
      ! 長い時間ステップのループを回すための処置
      !----------------------------------------------------------------

      xz_PotTempBl = xz_PotTempNl
      xza_MixRtBl  = xza_MixRtNl
      xz_ExnerBl   = xz_ExnerNl
      pz_VelXBl    = pz_VelXNl
      xr_VelZBl    = xr_VelZNl
      xz_KmBl      = xz_KmNl
      xz_KhBl      = xz_KhNl
      xz_DensCloudBl = xz_DensCloudNl
      xz_SatRatioBl = xz_SatRatioNl

      
      xz_PotTempNl = xz_PotTempAl
      xza_MixRtNl  = xza_MixRtAl
      xz_ExnerNl   = xz_ExnerAl
      pz_VelXNl    = pz_VelXAl
      xr_VelZNl    = xr_VelZAl
      xz_KmNl      = xz_KmAl
      xz_KhNl      = xz_KhAl
      xz_DensCloudNl = xz_DensCloudAl
      xz_SatRatioNl = xz_SatRatioAl

!      xz_PotTempSum = xz_PotTempNl + xz_PotTempBasicZ    
      xz_ExnerSum = xz_ExnerNl + xz_ExnerBasicZ
      xz_TempSum = xz_PotTempSum * xz_ExnerSum
        ! 2008/08/16(山下 達也)
        ! 基本場と擾乱場の和を出力する為に追加

   ! 積分量の計算

    ! 各格子点の質量密度(擾乱成分の寄与)
     xz_MassDens =                                                    &
      & (xz_ExnerBasicZ + xz_ExnerNl)**( CvDry/GasRDry )              &
      & / (xz_PotTempBasicZ + xz_PotTempNl)                           &
      & - (xz_ExnerBasicZ )**( CvDry/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))

    ! 温位の 2 次移流項を残すことによって生じてしまう
    ! 全質量への寄与
     xz_MassTend = - xz_DensBasicZ * xz_TendPotTempNl / xz_PotTempBasicZ

     MassTend = DelX * DelZ                                     &
      & * sum( xz_MassTend(RegXMin:RegXMax,RegZMin:RegZMax) )

    ! 各格子点の運動エネルギー(密度を基本場の密度で評価)
     xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry      &
      & * (xz_ExnerBasicZ)**( CvDry/GasRDry )                         &
      & / (xz_PotTempBasicZ)                                          &
      & * ( xz_avr_pz(pz_VelXNl)**2.0d0 + xz_avr_xr(xr_VelZNl)**2.0d0 &
      &     )

    ! 領域全体の運動エネルギー
     KineticEnergyTotal =                                             &
      & sum( xz_KineticEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )

    ! 領域全体の雲の質量
     CloudMassTotal = DelX * DelZ                                     &
      & * sum( xz_DensCloudNl(RegXMin:RegXMax,RegZMin:RegZMax) )

    ! 鉛直座標を 2 次元配列化
      do k = RegZMin, RegZMax
        xz_Z(:,k) = s_Z(k)
      end do

    ! 各格子点の位置エネルギー
     xz_PotentialEnergy = - xz_DensBasicZ * xz_PotTempNl * 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_ExnerNl

    ! 領域全体の弾性エネルギー(1 次量)
     ElasticEnergyFOTotal = DelX * DelZ                               &
      & * sum( xz_ElasticEnergyFO(RegXMin:RegXMax,RegZMin:RegZMax) )

    ! 各格子点の弾性エネルギー(2 次量)
     xz_ElasticEnergySO = 0.5d0 * xz_DensBasicZ *                     &
      & (CpDry * xz_PotTempBasicZ * xz_ExnerNl / xz_VelSoundBasicZ )**2.0d0

    ! 領域全体の弾性エネルギー(2 次量)
     ElasticEnergySOTotal = DelX * DelZ                               &
      & * sum( xz_ElasticEnergySO(RegXMin:RegXMax,RegZMin:RegZMax) )

    end do

   !----------------------------------------------------------------
    ! ファイル出力
    !----------------------------------------------------------------
    ! CFL 条件のチェック
    call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )
    
    call CFLCheckTimeLongVelX( pz_VelXNl )
    call CFLCheckTimeLongVelZ( xr_VelZNl )
    
    !ヒストリファイル出力
    call HistoryFile_Output(  &
      &  Time,                &
      &  xz_PotTempNl,        &
      &  xz_PotTempSum,       &
      &  xz_TempSum,          &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xza_MixRtNl,         &
      &  xz_KmNl,             &
      &  xz_KhNl,             &
      &  xz_DensCloudNl,      &
      &  xz_SatRatioNl,       &
      &  MassTotal,           &
      &  MassTend,            &
      &  KineticEnergyTotal,  &
      &  CloudMassTotal,      &
      &  PotentialEnergyTotal, &
      &  ElasticEnergyFOTotal, &
      &  ElasticEnergySOTotal  &
      & )
    
    !積算値のクリア
    call StorePotTempClean
    call StoreDensCloudClean
    call StoreMixRtClean
    call StoreBuoyClean
    call StoreStabClean
    
  end do
  
  !----------------------------------------------------------------
  ! 出力ファイルのクローズ
  !----------------------------------------------------------------
  call HistoryFile_Close


  !----------------------------------------------------------------
  ! リスタートファイルの作成
  !----------------------------------------------------------------
  call ReStartFile_Open( )
  call ReStartFile_OutPut(                            &
    & Time - DelTimeLong,                             &
    & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
    & xza_MixRtBl,  xz_KmBl,    xz_KhBl,              &
    & xz_DensCloudBl, xz_SatRatioBl                   )
  call ReStartFile_OutPut(                            &
    & Time,                                           &
    & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
    & xza_MixRtNl,  xz_KmNl,    xz_KhNl,              &
    & xz_DensCloudNl, xz_SatRatioNl                   )
  call ReStartFile_Close
  
contains

  subroutine ArareAlloc
    !
    !初期化として, 配列を定義し, 値としてゼロを代入する.
    !

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate(                                                   &
      & pz_VelXBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xr_VelZBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xz_ExnerBl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerSum(DimXMin:DimXMax, DimZMin:DimZMax),         &
!
      & xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_PotTempSum(DimXMin:DimXMax, DimZMin:DimZMax),       &
      & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNs(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAs(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_TempSum(DimXMin:DimXMax, DimZMin:DimZMax),          &
!
      & xz_KmBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xz_KhBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xza_MixRtWork(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
      & xza_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
!
      & pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_FluxDensCloud(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_TurbDensCloud(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_NumDensCloud(DimXMin:DimXMax, DimZMin:DimZMax),       &
      & xz_TendDensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax),    &
      & xz_TendPotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),      &
!
      & xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),  &
!
      & DelTimeLFrog(NstepLong), NStepEular(NStepLong),          &
!
      & xz_MassCondNs(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_MassCondNl(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_DensCloudBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_DensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_DensCloudAl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_DensCloudNs(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_DensCloudAs(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_DensCloudWorkAs(DimXMin:DimXMax, DimZMin:DimZMax),    &
      & xz_SatRatioBl(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_SatRatioNl(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_SatRatioAl(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_SatRatioNs(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_SatRatioAs(DimXMin:DimXMax, DimZMin:DimZMax),         &
      & xz_QCond(DimXMin:DimXMax, DimZMin:DimZMax),              &
      & xz_QRad(DimXMin:DimXMax, DimZMin:DimZMax),               &
      & xz_QDis(DimXMin:DimXMax, DimZMin:DimZMax),               &
      & worknum(DimXMin:DimXMax, DimZMin:DimZMax),               &
      & xz_MassDens(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xz_MassTend(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)                   )

    pz_VelXNs = 0.0d0;     pz_VelXAs = 0.0d0    
    xr_VelZNs = 0.0d0;     xr_VelZAs = 0.0d0
    xz_ExnerNs = 0.0d0;    xz_ExnerAs = 0.0d0
    xz_PotTempNs = 0.0d0;     xz_PotTempAs = 0.0d0

    pz_VelXBl = 0.0d0;     pz_VelXNl = 0.0d0;     pz_VelXAl = 0.0d0
    xr_VelZBl = 0.0d0;     xr_VelZNl = 0.0d0;     xr_VelZAl = 0.0d0
    xz_ExnerBl = 0.0d0;    xz_ExnerNl = 0.0d0;    xz_ExnerAl = 0.0d0
    xz_KmBl = 0.0d0;       xz_KmNl = 0.0d0;       xz_KmAl = 0.0d0
    xz_KhBl = 0.0d0;       xz_KhNl = 0.0d0;       xz_KhAl = 0.0d0
    xz_PotTempBl = 0.0d0;  xz_PotTempNl = 0.0d0;  xz_PotTempAl = 0.0d0
    xza_MixRtBl = 0.0d0;    xza_MixRtNl = 0.0d0;    xza_MixRtAl = 0.0d0

    pz_AccelVelXNl = 0.0d0
    xr_AccelVelZNl = 0.0d0

    xz_FluxDensCloud = 0.0d0
    xz_TurbDensCloud = 0.0d0
    xz_NumDensCloud = 0.0d0
    xz_TendDensCloudNl = 0.0d0
    xz_TendPotTempNl = 0.0d0

    xza_DelMixRt = 0.0d0

    xz_MassCondNs = 0.0d0;       xz_MassCondNl = 0.0d0
    xz_DensCloudBl = 0.0d0;      xz_DensCloudNl = 0.0d0
    xz_DensCloudAl = 0.0d0
    xz_DensCloudNs = 0.0d0;      xz_DensCloudAs = 0.0d0
    xz_DensCloudWorkAs = 0.0d0
    xz_SatRatioBl = 0.0d0;       xz_SatRatioNl = 0.0d0
    xz_SatRatioAl = 0.0d0
    xz_SatRatioNs = 0.0d0;       xz_SatRatioAs = 0.0d0
    xz_QCond = 0.0d0;            xz_QRad = 0.0d0
    xz_QDis = 0.0d0            

    DelTimeEular = 0.0d0
    DelTimeLFrog = 0.0d0
    NStepEular = 0.0d0
    NstepLFrog = 0.0d0
    
    worknum = 0.0d0
    xz_MassDens = 0.0d0
    xz_MassTend = 0.0d0
    xz_KineticEnergy = 0.0d0
    xz_PotentialEnergy = 0.0d0
    xz_ElasticEnergyFO = 0.0d0
    xz_ElasticEnergySO = 0.0d0
    xz_Z = 0.0d0

  end subroutine ArareAlloc
  
  
end program arare
