!= Program ArareMPI
!
! Authors::   SUGIYAMA Ko-ichiro (杉山耕一朗), ODAKA Masatsugu (小高正嗣)
! Version::   $Id: ararempi.f90,v 1.13 2011-03-01 03:04:02 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
! 非静力学モデル deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 方程式系は準圧縮系.
!
!== Future Plans
!
!

program ararempi
  !
  !非静力学モデル 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

  !  MPI の初期化
  use mpiset,        only : mpii_init, mpii_end, myrank

  !  入出力ファイル名管理モジュール
  use filesetmpi,    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

  !  基本場設定モジュール
  use basicset,      only : basicset_init, &
    &                       xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, &
    &                       xz_VelSoundBasicZ

  !  積算値管理モジュール
  use StorePotTemp,      only : StorePotTemp_init, StorePotTempClean, &
    &                           StorePotTempCond
  use StoreMixRt,        only : StoreMixRt_init, StoreMixRtClean, &
    &                           StoreMixRtCond, StoreMixRtFill1, StoreMixRtFill2
  use StoreBuoy,         only : StoreBuoy_init, StoreBuoyClean
  use StoreStab,         only : StoreStab_init, StoreStabClean
  use StoreMom,          only : StoreMom_Init, StoreMomClean

  !  湿潤ルーチン設定モジュール
  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 boundarympi,    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 RestartFileIOmpi, only : ReStartFile_Open, ReStartFile_OutPut, &
    &                       ReStartFile_Close, ReStartFile_Get

  !  ヒストリファイル入出力モジュール
  use HistoryFileIOmpi, 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 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,   xz_TurbKm
  
  !  放射強制計算用モジュール
  use Radiation,    only : Radiation_init,  &
    &                      xz_RadHeatConst  
  
  !  地表フラックス計算用モジュール
!  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, 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_PotTempWork(:,:)
  real(8), allocatable :: xz_PotTempAllBl(:,:)
  real(8), allocatable :: xz_PotTempAllNl(:,:)
  real(8), allocatable :: xz_PotTempBl(:,:)
  real(8), allocatable :: xz_PotTempNl(:,:)
  real(8), allocatable :: xz_PotTempAl(:,:)
  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_MixRtAllBl(:,:,:)
  real(8), allocatable :: xza_MixRtAllNl(:,:,:)
  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)              :: Time
  real(8)              :: ReStartTime(2)
  real(8), allocatable :: DelTimeLFrog(:)
  real(8)              :: DelTimeEular
  integer, allocatable :: NStepEular(:)
  integer              :: NStepLFrog

  real(8), allocatable :: Km1(:,:), Km2(:,:), Km3(:,:), Km4(:,:), Km5(:,:), Km6(:,:), Km7(:,:)

  integer :: & 
    & t,     & 
    & tau,   & ! do ループ変数 ; do loop variable  
    & t1,    & ! do ループ変数 ; do loop variable  
    & t2       ! do ループ変数 ; do loop variable 

  logical :: CloudFlag = .true.
  logical :: TurbFlag  = .true.


  !コマンドライン引数の解釈
  !
  call mpii_init

  ! 初期化手続き Initialize procedure 
  !
  call MainInit

  !----------------------------------------------------------------------
  ! 数値積分
  !----------------------------------------------------------------------
!  do t1 = 1, NstepLFrog / NstepDisp
!    do t2 = 1, NstepDisp
  do t = 1, NstepLFrog   

      !時刻の設定
      Time = Time + DelTimeLong
!      t = (t1 - 1) * NstepDisp + t2
      
      !----------------------------------------------------------------
      ! 渦粘性係数, 渦拡散係数を求める.
      !----------------------------------------------------------------
      if (TurbFlag) then 
      xz_KmAl =                                                       &
        & xz_KmBl                                                     &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)               &
        &     + xz_NumDiffKm(xz_KmBl)                                 &
        &     + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) &
        &     + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)             &   
        &     + xz_TurbKm(xz_KmBl)                                    &
        &     + xz_DispKm(xz_KmBl)                                    &
        &    )

      Km1 = xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)
      call StoreDKmDtAdv(Km1)
      Km2 = xz_NumDiffKm(xz_KmBl)
      call StoreDKmDtDiff(Km2)
      Km3 = xz_BuoyKm(xz_PotTempBl)
      call StoreDKmDtBuoyT(Km3)
      Km4 = xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) 
      call StoreDKmDtBuoyM(Km4 - km3)
      Km5 = xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)
      call StoreDKmDtShear(Km5)
      Km6 = xz_TurbKm(xz_KmBl)
      call StoreDKmDtTurb(Km6)
      Km7 = xz_DispKm(xz_KmBl)
      call StoreDKmDtDisp(Km7)

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

      else
        xz_KmAl = 0.0d0
        xz_KmNl = 0.0d0
      end if      

      !境界条件
      call BoundaryXCyc_xz( xz_KmAl )
      call BoundaryZSym_xz( xz_KmAl )
      
      !渦拡散定数を決める
      xz_KhAl = 3.0d0 * xz_KmAl
      
      !----------------------------------------------------------------
      ! 温位の移流計算.
      !----------------------------------------------------------------    
      !時間積分
      xz_PotTempAllNl = xz_PotTempNl + xz_PotTempBasicZ
      xz_PotTempAllBl = xz_PotTempBl + xz_PotTempBasicZ
      xz_PotTempAl =                                                  &
        &   xz_PotTempBl                                              &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + xz_AdvScalar( xz_PotTempAllNl,  pz_VelXNl, xr_VelZNl) &
        &     + xz_TurbScalar(xz_PotTempAllBl,  xz_KhBl)              &
        &     + xz_NumDiffScalar(xz_PotTempBl)                        &
        &     + xz_DispHeat( xz_KmBl )                                &
        &     + xz_RadHeatConst( xz_ExnerBl )                         &
        &     + xz_HeatFluxDiff( xz_PotTempNl )                       &
!!      &     + xz_HeatFluxBulk( xz_PotTempNl )                     &
!!      &     + xz_NewtonCool( xz_PotTempBl )                       &
        &      )

      !境界条件
      call BoundaryXCyc_xz( xz_PotTempAl )
      call BoundaryZSym_xz( xz_PotTempAl )
      
      !----------------------------------------------------------------
      ! 凝縮成分の混合比の移流計算.
      !----------------------------------------------------------------
      xza_MixRtAllNl = xza_MixRtNl + xza_MixRtBasicZ
      xza_MixRtAllBl = xza_MixRtBl + xza_MixRtBasicZ
      xza_MixRtAl =                                                & 
        &   xza_MixRtBl                                            &
        & + DelTimeLFrog(t)                                        &
        &   * (                                                    &
        &    + xza_AdvScalar(xza_MixRtAllNl,  pz_VelXNl, xr_VelZNl)&
        &    + xza_TurbScalar(xza_MixRtAllBl, xz_KhBl)             &
        &    + xza_NumDiffScalar(xza_MixRtBl)                      &
        &    + xza_FallRain(xza_MixRtBl)                           &
        &    + xza_MixRtFluxDiff(xza_MixRtNl)                      &
        &   )

      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !移流によって負になった部分を埋める
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = xza_FillNegative_xza( xza_MixRtWork ) 
      
      !埋めた/削った量を保管
      call StoreMixRtFill1( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )

      !境界条件
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )
      
!      ! Negative values mixing ratios are corrected. 
!      xza_MixRtWork = xza_MixRtAl
!      xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
!      call StoreMixRtFill2( &
!        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
!        & )

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

      !-------------------------------------------------------------
      ! 暖かい雨のパラメタリゼーション.
      !   雲<-->雨 の変換を行う
      !-------------------------------------------------------------
      !雲から雨への変換分を追加する. 
      !  xza_Cloud2Rain 関数は, 引数として時間刻みを取ることで, 
      !  時間積分値を出力する. 

      !これまでの値を作業配列に保管
      xza_MixRtWork = xza_MixRtAl
      
      !雨への変化量を計算
      xza_MixRtAl   = xza_MixRtWork &
        &             + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog(t) )
      
      !雲から雨への変換量を保管
      call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) 
      
      !境界条件
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! 湿潤飽和調節法
      !   水蒸気<-->雲の変換を行う.
      !-------------------------------------------------------------
      !これまでの値を作業配列に保管
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork  = xza_MixRtAl
      
      !湿潤調節法を適用
      call MoistAdjustSvapPress( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) )
      call MoistAdjustNH4SH( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) )
      
      !湿潤調節法による温位と混合比の変化量を保管
      call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) )
      call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) 

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

      !-------------------------------------------------------------
      ! 暖かい雨のパラメタリゼーション.
      !   蒸気<-->雨 の変換を行う
      !-------------------------------------------------------------
      !雨から蒸気への変換に伴う混合比の変化を計算
      !  xza_Rain2Gas 関数は, 引数として時間刻みを取ることで, 
      !  時間積分値を出力する. 
     
      !これまでの値を作業配列に保管
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork = xza_MixRtAl
      
      !雨から蒸気への混合比変化を求める
      !  温位の計算において, 混合比変化が必要となるため, 
      !  混合比変化を 1 つの配列として用意する.
      xza_DelMixRt = 0.0d0
      xza_DelMixRt =                                                    &
        & (                                                             &
        &   + xza_Rain2Gas(                                             &
        &        xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
        &       )                                                       &
        &   + xza_Rain2GasNH4SH(                                        &
        &        xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
        &       )                                                       &
        &  )    
     
      !温位の計算. 雨から蒸気への変換に伴う潜熱・反応熱を追加.
      xz_PotTempAl =                                                                      &
        & xz_PotTempWork                                                                  &
        & + (                                                                             &
        &      + xz_Rain2GasHeat(xz_PotTempAl, xz_ExnerNl, xza_DelMixRt, DelTimeLFrog(t)) &
        &      + xz_Rain2GasHeatNH4SH(xz_ExnerNl, xza_DelMixRt, DelTimeLFrog(t))          &
        &    )
     
      !混合比の計算. 雨から蒸気への変換分を追加
      xza_MixRtAl   = xza_MixRtWork + xza_DelMixRt
      
      !雨から蒸気の値の保管
      call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) ) 
      call StoreMixRtCond( xza_DelMixRt / DelTimeLFrog(t) ) 
     
      !境界条件
      call BoundaryXCyc_xz( xz_PotTempAl )
      call BoundaryZSym_xz( xz_PotTempAl )
      call BoundaryXCyc_xza( xza_MixRtAl )
      call BoundaryZSym_xza( xza_MixRtAl )

      !-------------------------------------------------------------
      ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力
      !-------------------------------------------------------------
      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

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

        !-------------------------------------------------------------
        ! 速度 u の計算
        !-------------------------------------------------------------
        pz_VelXAs = &
          & pz_VelXNs                                          &
          &  + DelTimeEular                                    &
          &    * (                                             &
          &     - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)  &
          &     + pz_AccelVelXNl                               &
          &     )
        
        !境界条件
        call BoundaryXCyc_pz( pz_VelXAs )
        call BoundaryZSym_pz( pz_VelXAs )
        
        !-------------------------------------------------------------
        ! エクスナー関数の計算
        !-------------------------------------------------------------
        xz_ExnerAs = xz_Exner( &
          & xr_AccelVelZNl,    &
          & pz_VelXNs,         &
          & pz_VelXAs,         &
          & xr_VelZNs,         &
          & xz_ExnerNs)
        
        !境界条件
        call BoundaryXCyc_xz( xz_ExnerAs )
        call BoundaryZSym_xz( xz_ExnerAs )
        
        !-------------------------------------------------------------
        ! 速度 w の計算
        !-------------------------------------------------------------
        xr_VelZAs =                                                    &
          & xr_VelZNs                                                  &
          &  + DelTimeEular                                            &
          &    * (                                                     &
          &     - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
          &     + xr_AccelVelZNl                                       &
          &     )
        
        !境界条件
        call BoundaryXCyc_xr( xr_VelZAs )
        call BoundaryZAntiSym_xr( xr_VelZAs )

        !-------------------------------------------------------------
        ! 短い時間ステップのループを回すための処置
        !-------------------------------------------------------------
        xz_ExnerNs  = xz_ExnerAs
        pz_VelXNs   = pz_VelXAs
        xr_VelZNs   = xr_VelZAs
        
      end do Eular

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

      !-------------------------------------------------------------
      ! アッセリンの時間フィルタ.  Nl の値をフィルタリング
      !-------------------------------------------------------------
      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_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl)
    
      !-------------------------------------------------------------
      ! スポンジ層
      !-------------------------------------------------------------
      call DampSponge_pz( pz_VelXAl,    pz_VelXBl,    DelTimeLFrog(t) )
      call DampSponge_xr( xr_VelZAl,    xr_VelZBl,    DelTimeLFrog(t) )
      call DampSponge_xz( xz_ExnerAl,   xz_ExnerBl,   DelTimeLFrog(t) )
      call DampSponge_xz( xz_PotTempAl, xz_PotTempBl, DelTimeLFrog(t) )
      call DampSponge_xz( xz_KmAl,      xz_KmBl,      DelTimeLFrog(t) )

      !
      ! Negative values mixing ratios are corrected. 
      !
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
      call StoreMixRtFill2( &
        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
        & )

      !--------------------------------------------------------------
      ! 解析値
      !---------------------------------------------------------------
      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_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

    end do
    !----------------------------------------------------------------
    ! ファイル出力
    !----------------------------------------------------------------
    ! CFL 条件のチェック
!    if (myrank == 0) then 
!      call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )
!    end if
!    call CFLCheckTimeLongVelX( pz_VelXNl )
!    call CFLCheckTimeLongVelZ( xr_VelZNl )
    
    !ヒストリファイル出力
    call HistoryFile_Output(  &
      &  Time,                &
      &  xz_PotTempNl,        &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xza_MixRtNl,         &
      &  xz_KmNl,             &
      &  xz_KhNl              &
      & )

    !積算値のクリア
    call StorePotTempClean
    call StoreMixRtClean
    call StoreBuoyClean
    call StoreStabClean
    call StoreMomClean

    !----------------------------------------------------------------
    ! リスタートファイルの作成
    !----------------------------------------------------------------
    call ReStartFile_Open( t1 )
    call ReStartFile_OutPut(                            &
      & Time - DelTimeLong,                             &
      & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
      & xza_MixRtBl,  xz_KmBl,    xz_KhBl               )
    call ReStartFile_OutPut(                            &
      & Time,                                           &
      & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
      & xza_MixRtNl,  xz_KmNl,    xz_KhNl               )
    call ReStartFile_Close
    
  end do
  
  !----------------------------------------------------------------
  ! 出力ファイルのクローズ
  !----------------------------------------------------------------
  call HistoryFile_Close

  !----------------------------------------------------
  ! MPI END
  !----------------------------------------------------
  call mpii_end

contains

  subroutine MainInit

    implicit none

    !NAMELIST ファイル名の読み込み
    if (myrank == 0) then 
      call MessageNotify( "M", "main", "Namelist file is '%c'", c1=trim(cfgfile) )
    end if
    call argset_init(cfgfile)
    
    !デバッグ設定
    call debugset_init(cfgfile, myrank)
    
    !物質特性の初期化
    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 StorePotTemp_init( )
    call StoreMixRt_init( )
    call StoreBuoy_init( )  
    call StoreStab_init( )  
    call StoreMom_init( )  
    
    !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
    call ArareAlloc
    
    !予報変数の初期化
    !  ReStartFile が設定されている場合にはファイルを読み込み, 
    !  設定されていない場合には, デフォルトの基本場と擾乱場を作る. 
    if (trim(InitFile) /= '') then    
      
      !基本場, 擾乱場の初期値を netCDF ファイルから取得する.
      call ReStartFile_Get(            &
        & ReStartTime,                 &
        & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
        & xza_MixRtBl,  xz_KmBl,    xz_KhBl,              &
        & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, &
        & xza_MixRtNl,  xz_KmNl,    xz_KhNl               )
      
    else
      !デフォルト設定の基本場, 擾乱場を作成する. 
      call BasicEnv()  ! 水平一様なので single CPU 版と同じルーチンを使える. 
      call DisturbEnvMPI(cfgfile,                            &
        & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, &
        & xza_MixRtBl,  xz_KmBl,    xz_KhBl               )
      
      ! 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
      
    end if
    
    !----------------------------------------------------------------------
    ! パッケージ型モジュールの初期化
    !   デフォルトの値から変更する必要のあるルーチンのみ初期化
    !----------------------------------------------------------------------
    call Damping_Init( cfgfile )      !波の減衰係数の初期化
    !  call NumDiffusion_Init3( 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_ExnerNl,          &
        &  pz_VelXNl,           &
        &  xr_VelZNl,           &
        &  xza_MixRtNl,         &
        &  xz_KmNl,             &
        &  xz_KhNl            )
    end if
    
    !時刻 0 の場合には出力
    if ( Time == 0 ) then
      call HistoryFile_Output(  &
        &  Time,                &
        &  xz_PotTempNl,        &
        &  xz_ExnerNl,          &
        &  pz_VelXNl,           &
        &  xr_VelZNl,           &
        &  xza_MixRtNl,         &
        &  xz_KmNl,             &
        &  xz_KhNl            )
    end if
    
    !----------------------------------------------------------------------
    ! 設定のチェック
    !----------------------------------------------------------------------
    !CFL 条件のチェック
    call CFLCheckTimeShort( xz_VelSoundBasicZ )
    
  end subroutine MainInit

!--------------------------------------------------------------------------
  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_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax),          &
!
      & xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_PotTempAllBl(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & xz_PotTempAllNl(DimXMin:DimXMax, DimZMin:DimZMax),     &
      & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAl(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_MixRtAllBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),&
      & xza_MixRtAllNl(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),        & 
!
      & xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),  &
!
      & DelTimeLFrog(NstepLong), NStepEular(NStepLong)               ) 

    pz_VelXNs = 0.0d0;     pz_VelXAs = 0.0d0    
    xr_VelZNs = 0.0d0;     xr_VelZAs = 0.0d0
    xz_ExnerNs = 0.0d0;    xz_ExnerAs = 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

    xza_DelMixRt = 0.0d0

    DelTimeEular = 0.0d0
    DelTimeLFrog = 0.0d0; 
    NStepEular = 0.0d0
    NstepLFrog = 0.0d0
    
  end subroutine ArareAlloc
    
end program ararempi
