!= deepconv/arare 湿潤大気対流計算用主プログラム
!
!= deepconv/arare main program for moist atmospheric convection
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare_earth.f90,v 1.1 2011-02-28 12:33:49 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

program arare
  !
  ! 非静力学モデル deepconv/arare 湿潤大気対流計算用主プログラム.
  !

  ! モジュール引用 ; use statement 
  !

  ! gtool5 関連 
  ! gtool5 modules
  !
  use dc_types,      only: STRING
  use dc_string,     only: StoA
  use dc_message,    only: MessageNotify

  ! 初期設定モジュール
  ! Initialize module
  !
  use argset,        only: argset_init  
  use fileset,       only: fileset_init, InitFile
  use debugset,      only: debugset_init
  use timeset,       only: timeset_init, DelTimeLong, DelTimeShort,  &
    &                      NstepDisp, NstepLong, NstepShort
  use gridset,       only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
    &                      SpcNum
  use basicset,      only: basicset_init, xza_MixRtBasicZ, xz_DensBasicZ, &
    &                      xz_PotTempBasicZ, xz_VelSoundBasicZ

  ! 化学量計算モジュール 
  ! Chemical calculation modules
  !
  use ChemCalc,      only: ChemCalc_init
  use chemdata,      only: chemdata_init

  ! 力学過程計算用関数モジュール
  ! Dynamical processes module
  !
  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
  
  ! 乱流拡散計算用モジュール
  ! Turbulent diffusion module
  !
  use Turbulence,    only: Turbulence_Init, xz_TurbScalar, xza_TurbScalar,  &
    &                      pz_TurbVelX, xr_TurbVelZ, xz_ShearKm, xz_DispKm, &
    &                      xz_DispHeat, xz_TurbKm, xz_BuoyKm
  
  ! 境界からのフラックス計算用モジュール
  ! Surface flux module
  !
  ! 地球用
  use HeatFlux_N1994, only: &
    &                 xz_HeatFluxBulk, xza_MixRtFluxBulk, pz_MomFluxBulk
  ! 木星用
  !  use HeatFlux, only: &
  !     &                 xz_HeatFluxDiff,  xza_MixRtFluxDiff

  ! 放射強制計算用モジュール
  ! Radiative forceing module
  !
  use Radiation,     only: Radiation_init2, xz_RadHeatConst, xz_RadHeatVary, xz_NewtonCool

  ! 湿潤過程計算用モジュール
  ! Moist processes modules
  !
  use moistset,      only: moistset_init
  use MoistAdjust,   only: MoistAdjustSvapPress
  use WarmRainPrm,   only: WarmRainPrm_Init, xz_Rain2GasHeat,   &
    &                      xza_Rain2Gas, xza_Cloud2Rain, xza_FallRain
  use MoistBuoyancy, only: MoistBuoy_Init, xz_BuoyMoistKm, xr_BuoyMolWt, &
    &                      xr_BuoyDrag
  use fillnegative,  only: FillNegative_init, xza_FillNegative_xza

  ! 安定度の計算
  ! Static stability calculation module
  !
  use ECCM,          only: ECCM_init, ECCM_Stab

  ! 数値拡散/摩擦計算用モジュール
  ! Numerical diffussion /dumping module
  !
  use NumDiffusion4th,only: NumDiffusion_Init, xz_NumDiffScalar, &
    &                       xz_NumDiffKm, xza_NumDiffScalar,     &
    &                       pz_NumDiffVelX, xr_NumDiffVelZ
  use damping,        only: damping_init, DampSponge_xz, DampSponge_pz, &
    &                       DampSponge_xr

  ! 積算値管理モジュール
  ! Monitor variables setup modules
  !
  use StorePotTemp,  only: StorePotTemp_init, StorePotTempClean, &
    &                      StorePotTempCond
  use StoreMixRt,    only: StoreMixRt_init, StoreMixRtClean, StoreMixRtCond, &
    &                      StoreMixRtFill1, StoreMixRtFill2
  use StoreMom,      only: StoreMom_init, StoreMomClean
  use StoreBuoy,     only: StoreBuoy_init, StoreBuoyClean
  use StoreStab,     only: StoreStab_init, StoreStabClean
  use StoreKm,       only: StoreKm_Init, StoreKmClean, StoreDKmDtAdv, &
    &                      StoreDKmDtTurb, StoreDKmDtDisp, StoreDKmDtDiff,  &
    &                      StoreDKmDtShear, StoreDKmDtBuoyT, StoreDKmDtBuoyM

  ! ファイル入出力モジュール
  ! File I/O module
  !
  use RestartFileIO, only: ReStartFile_Open, ReStartFile_OutPut, &
    &                      ReStartFile_Close, ReStartFile_Get
  use HistoryFileIO, only: HistoryFile_Open, HistoryFile_OutPut, &
    &                      HistoryFile_Close

  ! 下請けモジュール
  ! Utility modules
  !
  use cflcheck,      only: CFLCheckTimeShort, CFLCheckTimeLongVelX, &
    &                      CFLCheckTimeLongVelZ
  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

  ! 変数宣言 Variables definition 
  !

  implicit none

  ! 内部変数
  ! Internal variables
  !
  character(80) :: cfgfile  ! NAMELIST ファイル名 NAMELIST file name
  real(8), allocatable :: pz_VelXBl(:,:) 
                            ! $ u(t-\Delta t) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXNl(:,:) 
                            ! $ u(t) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXAl(:,:) 
                            ! $ u (t+\Delta t)a $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXNs(:,:) 
                            ! $ u (\tau) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXAs(:,:) 
                            ! $ u (\tau +\Delta \tau) 水平風 Horizontal wind
  real(8), allocatable :: xr_VelZBl(:,:) 
                            ! $ w (t-\Delta t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZNl(:,:) 
                            ! $ w (t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZAl(:,:) 
                            ! $ w (t+\Delta t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZNs(:,:) 
                            ! $ w (\tau) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZAs(:,:) 
                            ! $ w (\tau +\Delta \tau) $ 鉛直風 Vertical wind
  real(8), allocatable :: xz_ExnerBl(:,:)
                            ! $ \pi (t-\Delta t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerNl(:,:)
                            ! $ \pi (t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerAl(:,:)
                            ! $ \pi (t+\Delta t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerNs(:,:)
                            ! $ \pi (\tau) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerAs(:,:)
                            ! $ \pi (\tau +\Delta \tau) $ 圧力関数 Exner function
  real(8), allocatable :: xz_PotTempBl(:,:) 
                            ! $ \theta (t-\Delta t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempNl(:,:) 
                            ! $ \theta (t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempAl(:,:) 
                            ! $ \theta (t+\Delta t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempWork(:,:)
                            ! 温位の作業配列 Work array for potential temp.
  real(8), allocatable :: xz_KmBl(:,:)
                            ! $ Km (t-\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KmNl(:,:)         
                            ! $ K_m (t) 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KmAl(:,:)
                            ! $ K_m (t+\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhBl(:,:)
                            ! $ K_h (t-\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhNl(:,:)
                            ! $ K_h (t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhAl(:,:)
                            ! $ K_h (t+\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff.
  real(8), allocatable :: xza_MixRtBl(:,:,:)
                            ! $ q (t-\Delta t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtNl(:,:,:) 
                            ! $ q (t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtAl(:,:,:)
                            ! $ q (t+\Delta t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtWork(:,:,:) 
                            ! 湿潤量の作業配列
                            ! Work array for mixing ratio.
  real(8), allocatable :: pz_AccelVelXNl(:,:)  
                            ! 気圧傾度力を除いた $u$ の変化率
                            ! Tendency of $u$ except for pressure gradient term
  real(8), allocatable :: xr_AccelVelZNl(:,:)  
                            ! 気圧傾度力を除いた $w$ の変化率
                            ! Tendency of $w$ except for pressure gradient term
  real(8), allocatable :: xza_DelMixRt(:,:,:)  
                            ! 湿潤量の混合比 $ q $ の増分
                            ! Mixing ratio variation.
  real(8) :: Time           ! 時刻 Time 
  real(8) :: ReStartTime(2) ! リスタートファイル出力時刻用配列
                            ! Output time array for restart file
  real(8), allocatable :: DelTimeLFrog(:)
                            ! リープフロッグスキーム用時間格子間隔
                            ! Time interval for Leap-frog scheme
  real(8) :: DelTimeEular   ! オイラースキーム用時間格子
                            ! Time interval for Eular scheme
  integer :: NStepLFrog     ! リープフロッグスキーム用時間ステップ数
                            ! The number of time step for Leap-frog scheme
  integer, allocatable :: NStepEular(:)
                            ! オイラースキーム用時間ステップ数
                            ! The number of time step for Eular scheme

  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.


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

  ! 時間積分 time integration 
  !
!  do t1 = 1, NstepLFrog / NstepDisp 
!    do t2 = 1, NstepDisp 

  do t = 1, NstepLFrog

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

      ! 渦拡散係数の移流拡散
      ! Advection and diffusion of turbulent diffusion coefficient.
      !
      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, 学位論文)参照
      !
      ! Upper and lower bound value are specified.
      !
      xz_KmAl = max( 0.0d0, min( xz_KmAl, 800.0d0 ) )

      else
        xz_KmAl = 0.0d0
        xz_KmNl = 0.0d0
      end if

      ! 境界条件 Boundary condition
      call BoundaryXCyc_xz( xz_KmAl )  ! (inout)
      call BoundaryZSym_xz( xz_KmAl )  ! (inout)

      ! スカラーに対する渦拡散係数の計算 
      ! Specify turbulent diffusion coefficient for scalar variables.
      !
      xz_KhAl = 3.0d0 * xz_KmAl
      
      ! 温位の移流拡散の計算 
      ! Advection and diffusion of potential temperature.
      !
      xz_PotTempAl =                                                  &
        &   xz_PotTempBl                                              &
        & + DelTimeLFrog(t)                                           &
        &   * (                                                       &
        &     + 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_DispHeat( xz_KmBl )                                &
!!       &     + xz_RadHeatConst( xz_ExnerBl )                         &
!        &     + xz_RadHeatVary( xz_ExnerBl )                          &
!!       &     + xz_HeatFluxDiff( xz_PotTempNl )                       &
!        &     + xz_HeatFluxBulk( xz_PotTempNl, pz_VelXNl)             &
!        &     + xz_NewtonCool( xz_PotTempBl )                         &
        &      )

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xz( xz_PotTempAl ) ! (inout)
      call BoundaryZSym_xz( xz_PotTempAl ) ! (inout)
      
      ! 凝縮成分混合比の移流拡散 
      ! Advection and diffusion of vapor, cloud and rain mixing ratios.
      !
      xza_MixRtAl =                                                 & 
        &   xza_MixRtBl                                             &
        & + DelTimeLFrog(t)                                         &
        &   * (                                                     &
        &    + xza_AdvScalar(xza_MixRtNl,     pz_VelXNl, xr_VelZNl) &
        &    + xza_AdvScalar(xza_MixRtBasicZ, pz_VelXNl, xr_VelZNl) & 
        &    + xza_TurbScalar(xza_MixRtBl,    xz_KhBl)              &
        &    + xza_TurbScalar(xza_MixRtBasicZ,xz_KhBl)              &
        &    + xza_NumDiffScalar(xza_MixRtBl)                       &
        &    + xza_FallRain(xza_MixRtNl)                            &
!        &    + xza_FallRain(xza_MixRtBl)                            &
!!      &    + xza_MixRtFluxDiff(xza_MixRtNl)                       &
!!      &    + xza_MixRtFluxBulk(xza_MixRtNl)                       &
!!      &    + xza_MixRtFluxBulk(xza_MixRtNl, pz_VelXNl, xz_ExnerNl)&
!        &    + xza_MixRtFluxBulk(xza_MixRtNl, pz_VelXNl)         &
        &   )

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)

      ! 移流によって負になった部分を埋める
      ! Negative values due to advection are corrected.
      !
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = xza_FillNegative_xza( xza_MixRtWork ) 
      
      ! 埋めた/削った量を保管
      ! Correction value is stored.
      !
      call StoreMixRtFill1( &
        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
        & )    
      
      ! 埋め切れなかった部分をゼロにする. Fill2 に保管
      ! 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)
        & )

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)


      if (CloudFlag) then 

      ! 暖かい雨のパラメタリゼーション.
      ! * 雲<-->雨 の変換を行う.
      !
      ! Warm rain parameterization.
      ! * Conversion from cloud to rain.

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xza_MixRtWork = xza_MixRtAl
      
      ! 雨への変化量を計算
      ! Conversion values are calculated.
      !
      xza_MixRtAl = &
        & xza_MixRtWork + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog(t) )
      
      ! 雲から雨への変換量を保管
      ! Conversion values are sotred.
      !
      call StoreMixRtCond( &
        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
        & )
      
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)
      
      ! 湿潤飽和調節
      ! * 蒸気<-->雲の変換を行う.
      !
      ! Moist adjustment.
      ! * Conversion from vapor to cloud.

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork  = xza_MixRtAl
      
      ! 湿潤調節法を適用
      ! Moist adjustment is applied.
      !
      call MoistAdjustSvapPress( &
        & xz_ExnerNl,            & ! (in)
        & xz_PotTempAl,          & ! (inout)
        & xza_MixRtAl            & ! (inout)
        & )
      
      ! 湿潤調節法による温位と混合比の変化量を保管
      ! Adjustment values are stored.
      !
      call StorePotTempCond( &
        & (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t)  & ! (in)
        & )
      call StoreMixRtCond( &
        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) & ! (in)
        & )
      
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)

      ! 暖かい雨のパラメタリゼーション.
      ! * 蒸気<-->雨 の変換を行う
      !
      ! Warm rain parameterization.
      ! * Conversion from rain to vapor.

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork = xza_MixRtAl
      
      ! 雨から蒸気への混合比変化を求める
      ! * 温位の計算において, 混合比変化が必要となるため, 
      !   混合比変化を 1 つの配列として用意する.
      !
      ! Conversion values are calculated.
      !
      xza_DelMixRt = 0.0d0
      xza_DelMixRt =                                                   &
        & (                                                            &
        &   + xza_Rain2Gas(                                            &
        &       xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
        &       )                                                      &
        &  )    
      
      ! 温位の計算. 雨から蒸気への変換に伴う潜熱・反応熱を追加.
      !
      !
      xz_PotTempAl =                                                    &   
        & xz_PotTempWork                                                &
        & + (                                                           &
        &     xz_Rain2GasHeat( xz_PotTempAl, xz_ExnerNl, xza_DelMixRt ) & 
        &    )
      
      ! 混合比の計算. 雨から蒸気への変換分を追加
      !
      !
      xza_MixRtAl   = xza_MixRtWork + xza_DelMixRt
      
      ! 雨から蒸気の値の保管
      !
      !
      call StorePotTempCond( &
        & (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) & ! (in)
        & ) 
      call StoreMixRtCond( &
        & xza_DelMixRt / DelTimeLFrog(t) & ! (in)
        & ) 
 
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xz( xz_PotTempAl ) ! (inout)
      call BoundaryZSym_xz( xz_PotTempAl ) ! (inout)
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)

      end if
      
      ! 速度の移流拡散.
      ! Advection and diffusion of velocity components.
      !
      pz_AccelVelXNl =                                    &
        & + pz_AdvVelX(pz_VelXNl, xr_VelZNl)              &
        & + pz_TurbVelX(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
        & + pz_NumDiffVelX(pz_VelXBl)                     !&
!        & + pz_MomFluxBulk(pz_VelXNl)
     
      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)                        


      ! 短い時間ステップの初期値作成.
      ! Initial values set up for time integration with short time step.
      !
      pz_VelXNs  = pz_VelXBl
      xr_VelZNs  = xr_VelZBl
      xz_ExnerNs = xz_ExnerBl
      
      ! 短い時間ステップの時間積分. オイラー法を利用.
      ! Time integration with short time step.
      !
      Eular: do tau = 1, NstepEular(t)
        
        ! 速度 u の計算.
        ! Time integration horizontal velocity (u).
        !
        pz_VelXAs = &
          & pz_VelXNs                                          &
          &  + DelTimeEular                                    &
          &    * (                                             &
          &     - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)  &
          &     + pz_AccelVelXNl                               &
          &     )
        
        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_pz( pz_VelXAs ) ! (inout)
        call BoundaryZSym_pz( pz_VelXAs ) ! (inout)
        
        ! エクスナー関数の計算.
        ! Time integration exner function.
        !
        xz_ExnerAs = xz_Exner( &
          & xr_AccelVelZNl,    &
          & pz_VelXNs,         &
          & pz_VelXAs,         &
          & xr_VelZNs,         &
          & xz_ExnerNs)
        
        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_xz( xz_ExnerAs ) ! (inout)
        call BoundaryZSym_xz( xz_ExnerAs ) ! (inout)
        
        ! 速度 w の計算
        ! Time integration vertical velocity.
        !
        xr_VelZAs =                                                    &
          & xr_VelZNs                                                  &
          &  + DelTimeEular                                            &
          &    * (                                                     &
          &     - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
          &     + xr_AccelVelZNl                                       &
          &     )
        
        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_xr( xr_VelZAs )     ! (iuout)
        call BoundaryZAntiSym_xr( xr_VelZAs ) ! (inout)
        
        ! 短い時間ステップのループを回すための処置
        ! Renew prognostic variables for next short time step integration.
        !
        xz_ExnerNs  = xz_ExnerAs
        pz_VelXNs   = pz_VelXAs
        xr_VelZNs   = xr_VelZAs
        
      end do Eular
      
      ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
      ! Renew prognostic variables for next long time step integration.
      !
      xz_ExnerAl  = xz_ExnerAs
      pz_VelXAl   = pz_VelXAs
      xr_VelZAl   = xr_VelZAs
      
      ! 時間フィルタ. 
      ! Time filter. 
      !
      call AsselinFilter_xz( &
        &  xz_ExnerAl,       & ! (in)
        &  xz_ExnerNl,       & ! (inout)
        &  xz_ExnerBl        & ! (in)
        & )         
      call AsselinFilter_pz( &
        & pz_VelXAl,         & ! (in)
        & pz_VelXNl,         & ! (inout)
        & pz_VelXBl          & ! (in)
        & )
      call AsselinFilter_xr( &
        & xr_VelZAl,         & ! (in)
        & xr_VelZNl,         & ! (inout)
        & xr_VelZBl          & ! (in)
        & )
      call AsselinFilter_xz( &
        & xz_PotTempAl,      & ! (in)
        & xz_PotTempNl,      & ! (inout)
        & xz_PotTempBl       & ! (in)
        & )
      call AsselinFilter_xz( &
        & xz_KmAl,           & ! (in)
        & xz_KmNl,           & ! (inout)
        & xz_KmBl            & ! (in)
        & )
      call AsselinFilter_xza( &
        & xza_MixRtAl,        & ! (in)
        & xza_MixRtNl,        & ! (inout)
        & xza_MixRtBl         & ! (in)
        & )
    
      ! スポンジ層.
      ! Numerical dumping.
      !
      call DampSponge_pz( &
        & pz_VelXAl,      & ! (inout)
        & pz_VelXBl,      & ! (inout)
        & DelTimeLFrog(t) & ! (in)
        & )
      call DampSponge_xr( &
        & xr_VelZAl,      & ! (inout)
        & xr_VelZBl,      & ! (inout)
        & DelTimeLFrog(t) & ! (in)
        & )

      ! 安定度の計算.
      ! Calculation static stability.
      !
      call ECCM_Stab(   & 
        & xz_PotTempAl, & ! (in)
        & xz_ExnerAl,   & ! (in)
        & xza_MixRtAl   & ! (in)
        & )
      
      ! 長い時間ステップのループを回すための処置.
      ! Renew prognostic variables for next long time step integration.
      !
      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

    if ( mod((t - 1), NstepDisp) == 0 ) then 
    
    ! ヒストリーファイルへの出力.
    ! Out put to history file.
    !
    call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )

    ! 移流に対する CFL 条件のチェック
    ! CFL condtion check for advection
    !
    call CFLCheckTimeLongVelX( &
      & pz_VelXNl              & ! (in)
      & )
    call CFLCheckTimeLongVelZ( &
      & xr_VelZNl              & ! (in)
      & )
    
    ! ヒストリファイル出力
    ! Out put to history file.
    !
    call HistoryFile_Output( &
      & Time - DelTimeLong,  & ! (in)
      & xz_PotTempBl,        & ! (in) 
      & xz_ExnerBl,          & ! (in)
      & pz_VelXBl,           & ! (in)
      & xr_VelZBl,           & ! (in)
      & xza_MixRtBl,         & ! (in)
      & xz_KmBl,             & ! (in)
      & xz_KhBl              & ! (in)
      & )
    
    ! 積算値のクリア.
    ! Clear monitor variables.
    !
    call StorePotTempClean
    call StoreMixRtClean
    call StoreMomClean   ! 不要？
    call StoreBuoyClean
    call StoreStabClean
    call StoreKmClean

    end if
    
  end do
  
  ! 出力ファイルのクローズ
  ! Close out put files.
  !
  call HistoryFile_Close

  ! リスタートファイルの作成
  ! Make restartfile.
  !
  call ReStartFile_Open( )
  call ReStartFile_OutPut( &
    & Time - DelTimeLong,  & ! (in)
    & xz_PotTempBl,        & ! (in)
    & xz_ExnerBl,          & ! (in)
    & pz_VelXBl,           & ! (in)
    & xr_VelZBl,           & ! (in)
    & xza_MixRtBl,         & ! (in)
    & xz_KmBl,             & ! (in)
    & xz_KhBl              & ! (in)
    & )
  call ReStartFile_OutPut( &
    & Time,                &
    & xz_PotTempNl,        & ! (in)
    & xz_ExnerNl,          & ! (in)
    & pz_VelXNl,           & ! (in)
    & xr_VelZNl,           & ! (in)
    & xza_MixRtNl,         & ! (in)
    & xz_KmNl,             & ! (in)
    & xz_KhNl              & ! (in)
    & )
  call ReStartFile_Close
  
contains
  !-----------------------------------------------------------------------
  subroutine MainInit
    ! NAMELIST ファイル名の読み込み
    ! Loading NAMELIST file.
    !
    call argset_init( &
      & cfgfile       & ! (out) 
      & )
    
    ! デバッグ設定の初期化
    ! Initialization of debug output control.
    !
    call debugset_init( &
      & cfgfile         & ! (in)
      & )
    
    ! 化学定数の初期化
    ! Initialization of chemical constatns.
    !
    call chemdata_init()

    ! 時間積分の初期化
    ! Initialization of time integration.
    !
    call timeset_init( &
      & cfgfile        & ! (in)
      & )
    
    ! 格子点情報の初期化
    ! Initialization of grid arrangement.
    !
    call gridset_init( &
      & cfgfile        & ! (in)
      & )
    
    ! 化学計算ルーチンの初期化
    ! Initialization of chemical routines.
    !
    call chemcalc_init()
    
    ! 基本場変数の初期化
    ! Initialization of basic state variables.
    !
    call basicset_init( &
      & cfgfile         & ! (in)
      & )
    
    ! 湿度分布の初期化
    ! Initialization of moist profile variables.
    ! NAMELISTにeccmの新規項目追加に伴い編集(2010年12月)
    !
    call eccm_init( &
      & cfgfile         & ! (in)
      & )

    ! I/O ファイル名の初期化
    ! Initialization of output file name. 
    !
    call fileset_init( &
      & cfgfile        & ! (in)
      & )
    
    ! 湿潤過程共有変数の初期化
    ! Initialization of common variables for moist process.
    !
    call moistset_init( )
    
    ! 積算値保管変数の初期化
    ! Initialization of monitor variables.
    !
    call StorePotTemp_init( )
    call StoreMixRt_init( )
    call StoreMom_init( )
    call StoreBuoy_init( )  
    call StoreStab_init( )  
    call StoreKm_init( )  
    
    ! 内部変数の初期化
    ! Initialization of internal variables.
    !
    call ArareAlloc

    ! 初期値の代入 
    ! * ReStartFile が設定されている場合にはファイルを読み込む. 
    !   設定されていない場合にはデフォルトの基本場と擾乱場を作る. 
    !
    ! Initial value set up.
    ! * Read restartfile if it is specified. If not, make default basic
    !   state and disturbance.
    !
    call MessageNotify( "M", "main", "Initial value setup." )

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

      call MessageNotify( &
        & "M", "main", "Restart file is %c", c1=trim(Initfile) )
      
      call ReStartFile_Get( &
        & ReStartTime,      & ! (out)
        & xz_PotTempBl,     & ! (out)
        & xz_ExnerBl,       & ! (out)
        & pz_VelXBl,        & ! (out)
        & xr_VelZBl,        & ! (out)
        & xza_MixRtBl,      & ! (out)
        & xz_KmBl,          & ! (out)
        & xz_KhBl,          & ! (out)
        & xz_PotTempNl,     & ! (out)
        & xz_ExnerNl,       & ! (out)
        & pz_VelXNl,        & ! (out)
        & xr_VelZNl,        & ! (out)
        & xza_MixRtNl,      & ! (out)
        & xz_KmNl,          & ! (out)
        & xz_KhNl           & ! (out)
        & )
    else
      
      call MessageNotify( &
        & "W", "main", "Restart file is not specified." )

      call BasicEnv( )
      call DisturbEnv(  &
        & cfgfile,      & ! (out)
        & xz_PotTempBl, & ! (out)
        & xz_ExnerBl,   & ! (out)
        & pz_VelXBl,    & ! (out)
        & xr_VelZBl,    & ! (out)
        & xza_MixRtBl,  & ! (out)
        & xz_KmBl,      & ! (out)
        & xz_KhBl       & ! (out)
        & )
   
      call BoundaryXCyc_pz( pz_VelXBl )     ! (inout)
      call BoundaryZSym_pz( pz_VelXBl )     ! (inout)
      
      call BoundaryXCyc_xr( xr_VelZBl )     ! (inout)
      call BoundaryZAntiSym_xr( xr_VelZBl ) ! (inout)
      
      call BoundaryXCyc_xz( xz_PotTempBl )  ! (inout)
      call BoundaryZSym_xz( xz_PotTempBl )  ! (inout)
      
      call BoundaryXCyc_xza( xza_MixRtBl )  ! (inout)
      call BoundaryZSym_xza( xza_MixRtBl )  ! (inout)
      
    
      ! 時刻 $ t $ の変数値の初期値の設定
      ! * 1 ループ目だけは $ t $ の値を $ t-\Delta t$ の値と同じにする. 
      !   1 ステップ目はオイラー法で解く必要があるが, 1 ステップ目と
      !   それ以外のステップを別々にコーディングしたくない為
      !
      ! Set up initial value of time = "t" variables.
      !  
      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

    ! 数値摩擦係数の初期化
    ! Initialization of numerical friction coefficient.
    !
    call Damping_Init( &
      & cfgfile        & ! (in)
      & )      
  
    ! 数値拡散項の初期化
    ! Initialization of numerical diffusion term.
    ! NAMELISTにnumdiffusionの新規項目追加に伴い編集(2010年12月)
    !
    call NumDiffusion_Init( &
      & cfgfile             & ! (in)
      & )          
    
    ! 乱流拡散項の初期化
    ! Initialization of turbulent diffusion term.
    !
    call Turbulence_Init( )            
    
    ! 暖かい雨のパラメタリゼーションの初期化 
    ! Initialization of warmrain parameterization.
    !
    call WarmRainPrm_Init( &
      & cfgfile            & ! (in)
      & )  
    
    ! 負の湿潤量の補填計算の初期化
    ! Initialization of negative moist value correction.
    !
    call FillNegative_Init( &
      & xza_MixRtBasicZ,    & ! (in)
      & xz_DensBasicZ       & ! (in)
      & )
    
    ! 放射強制の初期化
    !  Initialization of radiative forcing.
    !
    call Radiation_Init2( &
      & cfgfile           & ! (in)    
      & )
    
    ! 湿潤気塊の浮力計算の初期化
    ! Initialization of moist buoyancy calculation.
    !
    call MoistBuoy_Init( ) 
    
    ! 圧力計算用係数行列の初期化 
    ! Initialization of coefficient matrix for exner function calculation.
    !
    call xz_Exner_Init( )              
    
    ! 時刻とループ回数の初期化
    ! Initialization of time integration.
    !
    NstepLFrog   = NstepLong 
    NstepEular   = NstepShort 
    DelTimeLFrog = DelTimeLong * 2.0d0 
    DelTimeEular = DelTimeShort
    
    ! 計算開始時刻と時間格子間隔の初期化
    ! * ReStartFile が設定されている場合, ファイルから読み込んだ値を指定.
    ! * ReStartFile が設定されてない場合
    !   * 開始時刻は 0.0
    !   * 1 ステップ目の時間格子間隔を陽に指定
    !
    ! Setup restart time and time interval. 
    ! * Read restartfile if it is specified.
    ! * If not, 
    !   * "t" is set to be 0.
    !   * Time intervals for 1st step are specified explicitly.
    !
    if ( trim(InitFile) /= '') then    
      Time = ReStartTime(2)               
    else
      Time = 0.0d0                    
      NstepEular(1)   = NstepShort /2 
      DelTimeLFrog(1) = DelTimeLong   
    end if
    
    ! ヒストリーファイルへの出力
    ! Out put to history file.
    !
    call HistoryFile_Open( )
    
    if ( trim(InitFile) /= '') then    
      call HistoryFile_Output( &
        & ReStartTime(2),      & ! (in)
        & xz_PotTempNl,        & ! (in)
        & xz_ExnerNl,          & ! (in)
        & pz_VelXNl,           & ! (in)
        & xr_VelZNl,           & ! (in)
        & xza_MixRtNl,         & ! (in)
        & xz_KmNl,             & ! (in)
        & xz_KhNl              & ! (in)
        & )
    end if
    
    if ( Time == 0 ) then
      call HistoryFile_Output( & ! (in)
        & Time,                & ! (in)
        & xz_PotTempNl,        & ! (in)
        & xz_ExnerNl,          & ! (in)
        & pz_VelXNl,           & ! (in)
        & xr_VelZNl,           & ! (in)
        & xza_MixRtNl,         & ! (in)
        & xz_KmNl,             & ! (in)
        & xz_KhNl              & ! (in)
        & )
    end if
    
    ! 音波に対する CFL 条件のチェック
    ! CFL condtion check for sound wave.
    !
    call CFLCheckTimeShort( &
      & xz_VelSoundBasicZ   & ! (in)
      & )
    
  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_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_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

    allocate(                                            &
         & Km1(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km2(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km3(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km4(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km5(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km6(DimXMin:DimXMax, DimZMin:DimZMax),        &
         & Km7(DimXMin:DimXMax, DimZMin:DimZMax)         &
         &) 
    Km1 = 0.0d0
    Km2 = 0.0d0
    Km3 = 0.0d0
    Km4 = 0.0d0
    Km5 = 0.0d0
    Km6 = 0.0d0
    Km7 = 0.0d0
    
  end subroutine ArareAlloc
!-----------------------------------------------------------------------
end program arare
