!= Program Arare
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare_3d.f90,v 1.1 2007/08/09 15:05:37 odakker Exp $
! Tag Name::  $Name: arare4-20070810 $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! 非静力学モデル deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 方程式系は準圧縮系.
!  * 移流れテスト計算用.
!
!== Future Plans
!
!

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

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

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

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

  !-----    管理モジュール   -----
  !  化学量計算モジュール
  !  化学量計算モジュール
  use ChemCalc, only: ChemCalc_init
  use chemdata, only: chemdata_init

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

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

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

  !  格子点管理モジュール 
  use gridset_3d,    only : gridset_init, &
    &                       DimXMin, DimXMax, &
    &                       DimYMin, DimYMax, &
    &                       DimZMin, DimZMax, &
    &                       SpcNum

  !  基本場設定モジュール
  use basicset_3d,   only : basicset_init, &
    &                       xyz_DensBasicZ, xyz_PotTempBasicZ, &
    &                       xyz_VelSoundBasicZ

  !  積算値管理モジュール
  use storeset_3d,   only : storeset_init, storeclean, &
       &                   StoreCond

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


  !-----    下請けモジュール   -----
  !  数値摩擦計算モジュール 
  use damping_3d,    only : damping_init, & 
    &                       xyz_DampSponge, pyz_DampSponge, xyr_DampSponge, &
    &                       xqz_DampSponge

  !  時間積分フィルターモジュール
  use timefilter_3d,    only : AsselinFilter

  !  境界条件適用モジュール
  use xyz_bc_module, only : BoundaryXCyc_xyz, BoundaryYCyc_xyz, &
    &                       BoundaryZSym_xyz, BoundaryZCyc_xyz, &
    &                       BoundaryXCyc_pyz, BoundaryYCyc_pyz, &
    &                       BoundaryZSym_pyz, BoundaryZCyc_pyz, &
    &                       BoundaryXCyc_xqz, BoundaryYCyc_xqz, &
    &                       BoundaryZSym_xqz, BoundaryZCyc_xqz, &
    &                       BoundaryXCyc_xyr, BoundaryYCyc_xyr, &
    &                       BoundaryZSym_xyr, BoundaryZCyc_xyr, &
    &                       BoundaryZAsym_xyr

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

  !  負の湿潤量の補填計算モジュール
  !use fillnegative,  only : FillNegative_init, xza_FillNegative_xza


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

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

  !-----       力学過程        -----
  !  力学過程計算用関数モジュール
  use DynFunc_3d,       only : xyz_AdvScalar, xyz_AdvScalar2, &
    &                       pyz_AdvVelX, xqz_AdvVelY,      &
    &                       xyr_Buoy, xyr_AdvVelZ, pyz_GradPi, xqz_GradPi
  
  !  力学過程陰解法計算用関数モジュール
  use DynImpFunc_3d,    only : xyz_Exner_init, xyz_Exner, xyr_GradPi
  
  !-----       物理過程        -----
  !  数値拡散計算用モジュール
  use NumDiffusion_3d,  only : NumDiffusion_Init, xyz_NumDiffScalar, &
    &                       pyz_NumDiffVelX, xqz_NumDiffVelY, xyr_NumDiffVelZ 

  !  乱流拡散計算用モジュール
  !use Turbulence,   only : Turbulence_Init, &
  !  &                      xz_TurbScalar, xza_TurbScalar, pz_TurbVelX, &
  !  &                      xr_TurbVelZ  , xz_ShearKm    , xz_DispKm,   &
  !  &                      xz_DispHeat

  !  放射強制計算用モジュール
  !use Radiation_3d,    only : Radiation_init,  &
  !  &                      xyz_RadHeatConst, xyz_NewtonCool
  
  !  地表フラックス計算用モジュール
  !use HeatFlux,     only : xz_HeatFluxBulk, xz_MixRtFluxBulk
  !use HeatFlux,     only : xz_HeatFluxDiff, xza_MixRtFluxDiff

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

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

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

  !  断熱上昇気塊の温度減率計算用モジュール
!  use ECCM_3D,       only : eccm_init 


  !暗黙の型宣言禁止
  implicit none

  !内部変数
  character(80) :: cfgfile
  Real(DP), allocatable :: pyz_VelXBl(:,:,:)
  real(DP), allocatable :: pyz_VelXNl(:,:,:)
  real(DP), allocatable :: pyz_VelXAl(:,:,:)
  real(DP), allocatable :: pyz_VelXNs(:,:,:)
  real(DP), allocatable :: pyz_VelXAs(:,:,:)
  Real(DP), allocatable :: xqz_VelYBl(:,:,:)
  real(DP), allocatable :: xqz_VelYNl(:,:,:)
  real(DP), allocatable :: xqz_VelYAl(:,:,:)
  real(DP), allocatable :: xqz_VelYNs(:,:,:)
  real(DP), allocatable :: xqz_VelYAs(:,:,:)
  real(DP), allocatable :: xyr_VelZBl(:,:,:)
  real(DP), allocatable :: xyr_VelZNl(:,:,:)
  real(DP), allocatable :: xyr_VelZAl(:,:,:)
  real(DP), allocatable :: xyr_VelZNs(:,:,:)
  real(DP), allocatable :: xyr_VelZAs(:,:,:)
  real(DP), allocatable :: xyz_ExnerBl(:,:,:)
  real(DP), allocatable :: xyz_ExnerNl(:,:,:)
  real(DP), allocatable :: xyz_ExnerAl(:,:,:)
  real(DP), allocatable :: xyz_ExnerNs(:,:,:)
  real(DP), allocatable :: xyz_ExnerAs(:,:,:)
  real(DP), allocatable :: xyz_PotTempWork(:,:,:)
  real(DP), allocatable :: xyz_PotTempBl(:,:,:)
  real(DP), allocatable :: xyz_PotTempNl(:,:,:)
  real(DP), allocatable :: xyz_PotTempAl(:,:,:)
  real(DP), allocatable :: xyz_KmBl(:,:,:)
  real(DP), allocatable :: xyz_KmNl(:,:,:)
  real(DP), allocatable :: xyz_KmAl(:,:,:)
  real(DP), allocatable :: xyz_KhBl(:,:,:)
  real(DP), allocatable :: xyz_KhNl(:,:,:)
  real(DP), allocatable :: xyz_KhAl(:,:,:)

  real(DP), allocatable :: pyz_AccelVelXNl(:,:,:)
  real(DP), allocatable :: xqz_AccelVelYNl(:,:,:)
  real(DP), allocatable :: xyr_AccelVelZNl(:,:,:)

  real(DP)              :: Time
  real(DP)              :: ReStartTime(2)
  real(DP), allocatable :: DelTimeLFrog(:)
  real(DP)              :: DelTimeEuler
  integer, allocatable  :: NStepEuler(:)
  integer               :: NStepLFrog

  integer               :: t, tau, t1, t2   


  !----------------------------------------------------------------------
  ! 変数参照型モジュールの初期化
  !----------------------------------------------------------------------

  !コマンドライン引数の解釈
  !  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 moistset_init()

  !積算値を保管するためのモジュールの初期化
  !  NAMELIST から情報を得て, 基本場を設定する.
  call storeset_init( )
  
  !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
  call VariableAllocate
  
  !予報変数の初期化
  !  ReStartFile が設定されている場合にはファイルを読み込み, 
  !  設定されていない場合には, デフォルトの基本場と擾乱場を作る. 
  if (trim(InitFile) /= '') then    

    !基本場, 擾乱場の初期値を netCDF ファイルから取得する.
    call ReStartFile_Get(            &
      & ReStartTime,                 &
      & xyz_PotTempBl, xyz_ExnerBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, &
      & xyz_KmBl,      xyz_KhBl,              &
      & xyz_PotTempNl, xyz_ExnerNl, pyz_VelXNl, xqz_VelYBl, xyr_VelZNl, &
      & xyz_KmNl,    xyz_KhNl                         )
  else
    !デフォルト設定の基本場, 擾乱場を作成する. 
    call BasicEnv_3d()
    call DisturbEnv_3d(cfgfile,                             &
      & xyz_PotTempBl, xyz_ExnerBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, &
      & xyz_KmBl,    xyz_KhBl                 )

    call BoundaryXCyc_pyz( pyz_VelXBl )
    call BoundaryYCyc_pyz( pyz_VelXBl )
    call BoundaryZCyc_pyz( pyz_VelXBl )

    call BoundaryXCyc_xqz( xqz_VelYBl )
    call BoundaryYCyc_xqz( xqz_VelYBl )
    call BoundaryZCyc_xqz( xqz_VelYBl )

    call BoundaryXCyc_xyr( xyr_VelZBl )
    call BoundaryYCyc_xyr( xyr_VelZBl )
    call BoundaryZCyc_xyr( xyr_VelZBl )

    call BoundaryXCyc_xyz( xyz_PotTempBl )
    call BoundaryYCyc_xyz( xyz_PotTempBl )
    call BoundaryZCyc_xyz( xyz_PotTempBl )

    
    ! 1 ループ目だけは, bl と nl の値を同じにしておく. 
    ! 1 ステップ目はオイラー法で解く必要があるが, 
    ! 1 ステップ目とそれ以外のステップを別々にコーディングしたくない為
    xyz_PotTempNl = xyz_PotTempBl
    xyz_ExnerNl   = xyz_ExnerBl
    pyz_VelXNl    = pyz_VelXBl
    xqz_VelYNl    = xqz_VelYBl
    xyr_VelZNl    = xyr_VelZBl
    xyz_KmNl      = xyz_KmBl
    xyz_KhNl      = xyz_KhBl

  end if


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

     

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

  !----------------------------------------------------------------
  ! ファイル入出力
  !----------------------------------------------------------------
  !ファイルオープン

  call HistoryFile_Open( )

  if ( trim(InitFile) /= '') then    
    call HistoryFile_Output(  &
      &  ReStartTime(2),      &
      &  xyz_PotTempNl,        &
      &  xyz_ExnerNl,          &
      &  pyz_VelXNl,           &
      &  xqz_VelYNl,           &
      &  xyr_VelZNl,           &
      &  xyz_KmNl,             &
      &  xyz_KhNl            )
  end if

  !時刻 0 の場合には出力
  if ( Time == 0 ) then
    call HistoryFile_Output(  &
      &  Time,                &
      &  xyz_PotTempNl,        &
      &  xyz_ExnerNl,          &
      &  pyz_VelXNl,           &
      &  xqz_VelYNl,           &
      &  xyr_VelZNl,           &
      &  xyz_KmNl,             &
      &  xyz_KhNl            )
  end if
  
  !----------------------------------------------------------------------
  ! 設定のチェック
  !----------------------------------------------------------------------
  !CFL 条件のチェック
  call CFLCheckTimeShort( xyz_VelSoundBasicZ )

  !----------------------------------------------------------------------
  ! 数値積分
  !----------------------------------------------------------------------
  call MessageNotify( "M", "main", "Time Integration Start" )

  do t1 = 1, NstepLFrog / NstepDisp
    do t2 = 1, NstepDisp

      t = (t1 - 1) * NstepDisp + t2

      
      !時刻の設定
      Time = Time + DelTimeLong
      call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )
      
      !----------------------------------------------------------------
      ! 温位の移流計算.
      !----------------------------------------------------------------    
      !時間積分
      xyz_PotTempAl =                                              &
        &   xyz_PotTempBl                                          &
        & + DelTimeLFrog(t)                                        &
        &   * (                                                    &
        &     + xyz_AdvScalar( xyz_PotTempNl,                      &
        &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) &
!!      &     + xyz_AdvScalar( xyz_PotTempBasicZ,                  &
!!      &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) &
!!      &     + xyz_NumDiffScalar(xyz_PotTempBl)                   &
!!      &     + xyz_TurbScalar(xyz_PotTempBl,     xyz_KhBl)        &
!!      &     + xyz_TurbScalar(xyz_PotTempBasicZ, xyz_KhBl)        &
!!      &     + xyz_DispHeat( xyz_KmBl )                           &
!!      &     + xyz_RadHeatConst( xyz_ExnerBl )                    &
!!      &     + xyz_HeatFluxDiff( xyz_PotTempNl )                  &
!!      &     + xyz_HeatFluxBulk( xyz_PotTempNl )                  &
!!      &     + xyz_NewtonCool( xyz_PotTempBl )                    &
        &      )

      !境界条件
      call BoundaryXCyc_xyz( xyz_PotTempAl )
      call BoundaryYCyc_xyz( xyz_PotTempAl )
      call BoundaryZCyc_xyz( xyz_PotTempAl )
!      call BoundaryZSym_xyz( xyz_PotTempAl )
      
      !-------------------------------------------------------------
      ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力
      !-------------------------------------------------------------

      pyz_AccelVelXNl = 0.0d0
!      pyz_AccelVelXNl =                                     &
!        & + pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)  
!        & + pyz_NumDiffVelX(pyz_VelXBl)                      
!        & + pyz_TurbVelX(xyz_KmBl,   pyz_VelXBl, xyr_VelZBl) 

      xqz_AccelVelYNl = 0.0d0
!      xqz_AccelVelYNl =                                      &
!        & + xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)  
!        & + xqz_NumDiffVelY(xqz_VelYBl)                      
!        & + xqz_TurbVelX(xyz_KmBl,   pyz_VelXBl, xyr_VelZBl) 

      xyr_AccelVelZNl = 0.0d0
!      xyr_AccelVelZNl =                                       &
!        & + xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)   &
!        & + xyr_Buoy(xyz_PotTempNl)                              
!        & + xyr_NumDiffVelZ(xyr_VelZBl)                      &  
!        & + xyr_TurbVelZ(xyz_KmBl,   pyz_VelXBl, xyr_VelZBl) 

      !-------------------------------------------------------------
      ! 短い時間ステップの初期値作成
      !-------------------------------------------------------------
      pyz_VelXNs  = pyz_VelXBl
      xqz_VelYNs  = xqz_VelYBl
      xyr_VelZNs  = xyr_VelZBl
      xyz_ExnerNs = xyz_ExnerBl
    
      !-------------------------------------------------------------
      ! 短い時間ステップの積分. オイラー法を利用
      !-------------------------------------------------------------
      Euler: do tau = 1, NstepEuler(t)

        !-------------------------------------------------------------
        ! 速度 u の計算
        !-------------------------------------------------------------
        pyz_VelXAs = &
          & pyz_VelXNs                                          &
          &  + DelTimeEuler                                     &
          &    * (                                              &
          &     - pyz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs) &
          &     + pyz_AccelVelXNl                               &
          &     )

        !境界条件
        call BoundaryXCyc_pyz( pyz_VelXAs )
        call BoundaryYCyc_pyz( pyz_VelXAs )
!        call BoundaryZSym_pyz( pyz_VelXAs )
        call BoundaryZCyc_pyz( pyz_VelXAs )

        !-------------------------------------------------------------
        ! 速度 v の計算
        !-------------------------------------------------------------
        xqz_VelYAs = &
          & xqz_VelYNs                                          &
          &  + DelTimeEuler                                     &
          &    * (                                              &
          &     - xqz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs) &
          &     + xqz_AccelVelYNl                               &
          &     )

        !境界条件
        call BoundaryXCyc_xqz( xqz_VelYAs )
        call BoundaryYCyc_xqz( xqz_VelYAs )
!        call BoundaryZSym_xqz( xqz_VelYAs )
        call BoundaryZCyc_xqz( xqz_VelYAs )

        !-------------------------------------------------------------
        ! エクスナー関数の計算
        !-------------------------------------------------------------
        xyz_ExnerAs = xyz_Exner( &
          & xyr_AccelVelZNl,     &
          & pyz_VelXNs, pyz_VelXAs,      &
          & xqz_VelYNs, xqz_VelYAs,      &
          & xyr_VelZNs, xyz_ExnerNs)


        !境界条件
        call BoundaryXCyc_xyz( xyz_ExnerAs )
        call BoundaryYCyc_xyz( xyz_ExnerAs )
!        call BoundaryZSym_xyz( xyz_ExnerAs )
        call BoundaryZCyc_xyz( xyz_ExnerAs )

        !-------------------------------------------------------------
        ! 速度 w の計算
        !-------------------------------------------------------------
        xyr_VelZAs =                                                   &
          & xyr_VelZNs                                                 &
          &  + DelTimeEuler                                            &
          &    * (                                                     &
          &     - xyr_GradPi(xyz_ExnerAs,xyz_ExnerNs,                  &
          &                  pyz_VelXNs,xqz_VelYNs, xyr_VelZNs)        &
          &     + xyr_AccelVelZNl                                      &
          &     )

        !境界条件
        call BoundaryXCyc_xyr( xyr_VelZAs )
        call BoundaryYCyc_xyr( xyr_VelZAs )
        call BoundaryZCyc_xyr( xyr_VelZAs )
!        call BoundaryZAsym_xyr( xyr_VelZAs )
      
        !-------------------------------------------------------------
        ! 短い時間ステップのループを回すための処置
        !-------------------------------------------------------------
        xyz_ExnerNs  = xyz_ExnerAs
        pyz_VelXNs   = pyz_VelXAs
        xqz_VelYNs   = xqz_VelYAs
        xyr_VelZNs   = xyr_VelZAs

      end do Euler
    
      !----------------------------------------------------------------
      ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
      !----------------------------------------------------------------
      xyz_ExnerAl  = xyz_ExnerAs
      pyz_VelXAl   = pyz_VelXAs
      xqz_VelYAl   = xqz_VelYAs
      xyr_VelZAl   = xyr_VelZAs

      !-------------------------------------------------------------
      ! アッセリンの時間フィルタ.  Nl の値をフィルタリング
      !-------------------------------------------------------------
      call AsselinFilter(xyz_ExnerAl, xyz_ExnerNl, xyz_ExnerBl)
      call AsselinFilter(pyz_VelXAl,  pyz_VelXNl,  pyz_VelXBl )
      call AsselinFilter(xqz_VelYAl,  xqz_VelYNl,  xqz_VelYBl )
      call AsselinFilter(xyr_VelZAl,  xyr_VelZNl,  xyr_VelZBl )
      call AsselinFilter(xyz_PotTempAl, xyz_PotTempNl, xyz_PotTempBl)
!      call AsselinFilter(xyz_KmAl, xyz_KmNl, xyz_KmBl)

      !-------------------------------------------------------------
      ! スポンジ層
      !-------------------------------------------------------------
!      pyz_VelXAl = pyz_DampSponge( pyz_VelXAl, pyz_VelXBl, DelTimeLFrog(t) )
!      xqz_VelYAl = xqz_DampSponge( xqz_VelYAl, xqz_VelYBl, DelTimeLFrog(t) )
!      xyr_VelZAl = xyr_DampSponge( xyr_VelZAl, xyr_VelZBl, DelTimeLFrog(t) )


      !----------------------------------------------------------------
      ! 長い時間ステップのループを回すための処置
      !----------------------------------------------------------------
      xyz_PotTempBl = xyz_PotTempNl
      xyz_ExnerBl   = xyz_ExnerNl
      pyz_VelXBl    = pyz_VelXNl
      xqz_VelYBl    = xqz_VelYNl
      xyr_VelZBl    = xyr_VelZNl
      xyz_KmBl      = xyz_KmNl
      xyz_KhBl      = xyz_KhNl
    
      xyz_PotTempNl = xyz_PotTempAl
      xyz_ExnerNl   = xyz_ExnerAl
      pyz_VelXNl    = pyz_VelXAl
      xqz_VelYNl    = xqz_VelYAl
      xyr_VelZNl    = xyr_VelZAl
      xyz_KmNl      = xyz_KmAl
      xyz_KhNl      = xyz_KhAl

    end do

    !----------------------------------------------------------------
    ! ファイル出力
    !----------------------------------------------------------------
    
    call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )

    ! CFL 条件のチェック
    call CFLCheckTimeLongVelX( pyz_VelXNl )
    call CFLCheckTimeLongVelY( xqz_VelYNl )
    call CFLCheckTimeLongVelZ( xyr_VelZNl )
    
    !ヒストリファイル出力
    
    call HistoryFile_Output(  &
      &  Time,                &
      &  xyz_PotTempNl,       &
      &  xyz_ExnerNl,         &
      &  pyz_VelXNl,          &
      &  xqz_VelYNl,          &
      &  xyr_VelZNl,          &
      &  xyz_KmNl,            &
      &  xyz_KhNl             &
      & )

    !積算値のクリア
    call StoreClean
    
  end do

  !----------------------------------------------------------------
  ! 出力ファイルのクローズ
  !----------------------------------------------------------------
  call HistoryFile_Close


contains

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

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate(                                                   &
      & pyz_VelXBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xqz_VelYBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyr_VelZBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_ExnerBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_PotTempWork(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_KmBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KmNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KmAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_KhBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KhNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KhAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & pyz_AccelVelXNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_AccelVelYNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_AccelVelZNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),& 
      & DelTimeLFrog(NstepLong), NStepEuler(NStepLong)               ) 

    pyz_VelXNs  = 0.0d0;    pyz_VelXAs = 0.0d0    
    xqz_VelYNs  = 0.0d0;    xqz_VelYAs = 0.0d0    
    xyr_VelZNs  = 0.0d0;    xyr_VelZAs = 0.0d0
    xyz_ExnerNs = 0.0d0;    xyz_ExnerAs = 0.0d0

    pyz_VelXBl  = 0.0d0;    pyz_VelXNl  = 0.0d0;    pyz_VelXAl  = 0.0d0
    xqz_VelYBl  = 0.0d0;    xqz_VelYNl  = 0.0d0;    xqz_VelYAl  = 0.0d0
    xyr_VelZBl  = 0.0d0;    xyr_VelZNl  = 0.0d0;    xyr_VelZAl  = 0.0d0
    xyz_ExnerBl = 0.0d0;    xyz_ExnerNl = 0.0d0;    xyz_ExnerAl = 0.0d0
    xyz_KmBl    = 0.0d0;    xyz_KmNl    = 0.0d0;    xyz_KmAl    = 0.0d0
    xyz_KhBl    = 0.0d0;    xyz_KhNl    = 0.0d0;    xyz_KhAl    = 0.0d0
    xyz_PotTempBl = 0.0d0;  xyz_PotTempNl = 0.0d0;  xyz_PotTempAl = 0.0d0

    pyz_AccelVelXNl = 0.0d0
    xqz_AccelVelYNl = 0.0d0
    xyr_AccelVelZNl = 0.0d0

!    xyza_DelMixRt = 0.0d0

    DelTimeEuler = 0.0d0
    DelTimeLFrog = 0.0d0; 
    NStepEuler = 0.0d0
    NstepLFrog = 0.0d0

  end subroutine VariableAllocate
      
end program arare
